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ADDITION OF HIGHER ORDER PLATE AND SHELL ELEMENTS 
INTO NASTRAN COMPUTER PROGRAM 

By 

R. Narayanaswami 1 
ABSTRACT 

Two higher order plate elements, the linear strain triangular membrane 
element TRIM6 and the qumtic bending element TRPLT1, and a shallow shell 
element TRSHL, suitable for inclusion into the NASTRAN® (NASA STRu ctural 
ANalysis) program are described m this report. Additions to the NASTRAN 
Theoretical Manual [NASA SP-221 (03)], the NASTRAN Users' Manual [NASA SP-222 
(03)], the NASTRAN Programmers' Manual [NASA SP-223 (03)], and the NASTRAN 
Demonstration Problem Manual [NASA SP-224 (03) ] , for inclusion of these 
elements into the NASTRAN program are presented herein 

INTRODUCTION 

New higher order plate (nonconforming qumtic) and shell elements 
suitable for inclusion into the NASTRAN (NASA STRu ctural ANalysis) program 
have been developed by Narayanaswami (refs. 1 and 2) . The linear strain 
triangular membrane element developed by Argyris (ref. 3) has proved to be an 
accurate element for solving for membrane action m plates. Preliminary 
studies m the NASTRAN environment indicate that these elements are more 
efficient than the existing plate and shell elements in NASTRAN. These 
elements have now been added to the level 16.0 version of NASTRAN® and the 
complete computer listing for the addition has been delivered to the NASTRAN 
Systems Management Office, NASA Langley Research Center, Hampton, Virginia 
This report describes the theoretical formulations, information pertaining to 
their use, programming details and demonstration problems pertaining to the 

1 Research Associate, Old Dominion University Research Foundation, Norfolk, 
Virginia 23508 



three elements, TRIM6, TRPLT1 and TRSHL m update form suitable for incorporation 
into the respective NASTRAN manuals (refs. 4, 5, 6 and 7). 

DETAILS OF ELEMENT FORMULATIONS 

The Linear Strain Triangular Membrane 
Element TRIM6 

First proposed by Argyris (ref. 3), this element has six nodes, three at 
the comers and three at the midpoints of the sides. The element uses a 
quadratic displacement field. The thickness of the element as well as the 
temperature distribution within the element are permitted to have bilinear 
variation; the three constants of the bilinear equation for the same are 
evaluated by the respective user-specified values at the three comer nodes 
of the element. The FORTRAN subroutines for stiffness, mass, thermal load 
vector, and stress data recovery have been coded and tested out m stand-alone 
computer programs. The updates for incorporating the element into the NASTRAN 
program have been prepared and checked out m NASTRAN® Level 16.0 versions. 

The element is currently designed for use with the statics and normal modes 
rigid formats of NASTRAN. 

The Higher Order Triangular Bending 
Element TRPLT1 

This element was developed by Narayanaswami (refs. 1 and 7) as a 
modification of the high precision triangular plate bending element developed 
by Cowper et al. (ref. 8). The element has six nodes, three at the 
corners and three at the midpoints of the sides. A qumtic displacement field 
is chosen for the transverse displacement. Transverse shear flexibility is 
taken into account m the stiffness formulation. The thickness of the 
element is permitted to have bilinear variations, the three constants of the 
bilinear equation for the same are evaluated by the respective user-specified 
values at the three comer nodes of the element. The FORTRAN subroutines for 
stiffness, mass, thermal load vector and stress data recovery have been coded 
and tested out m stand-alone computer programs. The updates for incorporating 
the element into the NASTRAN program have been prepared and checked out m 
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NASTRAN'"' Level 16.0 version. The element is currently designed for use with 
the statics and normal modes rigid formats of NASTRAN. 

The Triangular Shallow Shell 
Element TRSHL 

This element was developed by Narayanaswami (ref. 2) . In the element 
coordinate system, the element has 30 degrees of freedom (d.o f. ), viz , the 
three translations u , v , w in the x , y , z directions and the 

2 rotations a and g about the xz and yz planes, at each of the 

3 corner nodes and 3 midside nodes of the triangle. The membrane behavior of 
this element is approximated by the TRIM6 element, the bending behavior is 
approximated by the TRPLT1 element and the membrane-bending coupling is 
approximated using shallow shell theory of Novozhilov (ref. 8) The element 
is currently designed for use with the statics, normal modes and buckling 
rigid formats of NASTRAN. 

ADDITIONS OR MODIFICATIONS TO NASTRAN MANUALS 
AND SOURCE CODE 

The updates to the NASTRAN Manuals for the addition of these elements 
are given m Appendixes A, B, C and D. The list of subroutines that are 
being modified or added to the NASTRAN Source Code is given in Appendix E. 

CONCLUDING REMARKS 

The addition of higher order plate elements, TRIM6 and TRPLT1, and the 
triangular shallow shell element TRSHL, into the NASTRAN program is completed 
These elements are added to the Level 16.0 version of NASTRAN® . The 
demonstration problems indicate the excellent accuracy of these elements for 
solving plate and shell problems. The availability of these elements m 
NASTRAN enchances the program’s capability m these areas. 
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APPENDIX A 

Updates to the NASTRAN Theoretical Manual 
for the addition of TRIM6, TRPLT1, 
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5.8 PLATES 


NASTRAN includes two different shapes of plate elements (triangular 
and quadrilateral) and two different stress systems (membrane and bending) 
which are, at present, uncoupled. There are m all a total of 13 different 
forms of plate elements as follows: 

1. TRMEM - A triangular element with finite inplane stiffness and zero 
bending stiffness. 

2. TRIM6 - A triangular element with finite inplane stiffness and zero 
bending stiffness. Uses quadratic polynomial representation for membrane 
displacements; bilinear variation m terms of the planar coordinates is 
permitted for the thickness as well as the temperature distribution of the 
element 

3. TRBSC - The basic unit from which the bending properties of the other 
plate elements are formed. In stand-alone form, it is used mainly as a 
research tool. 

4. TRPLT - A triangular element with zero inplane stiffness and finite 
bending stiffness. It si composed of three basic bending triangles that are 
coupled to form a Clough composite triangle; see section 5 8 3.3 

5. TRPLT1 - Higher order bending element - a triangular element with 
zero mplane stiffness and finite bending stiffness. Uses qumtic polynomial 
representation for transverse displacement, bilinear variation in terms of 
the planar coordinates of the element is permitted for the element thickness 

6. TRIA1 - A triangular element with both mplane and bending stiffness 
It is designed for sandwich plates m which different materials can be 
referenced for membrane, bending, and transverse shear properties. 

7 . TRIA2 - A triangular element with both mplane and bending stiffness 
that assumes a solid homogeneous cross section 

8. QDMEM - A quadrilateral membrane element consisting of four over- 
lapping TRMEM elements 

9 QDMEM1 - An isoparametric quadrilateral membrane element. 
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7 



10. QDMEM2 - A quadrilateral membrane element consisting of four 
nonoverlapping TRMEM elements. 

11. QDPLT - A quadrilateral bending element. It is composed of four 
basic bending triangles. 

12. QUAD1 - A quadrilateral element with both inplane and bending 
stiffness , similar to TRIA1. 

13. QUAD2 - A quadrilateral element similar to TRIA2. 

Anisotropic material properties may be employed in all plate elements. TRMEM 
and TRBSC are the basic plate elements from which TRPLT, TRIA1, TRIA2 , 

QDMEM, QDMEM2, QDPLT, QUAD1, and QUAD2 elements are formed. The stiffness 
matrices of plate elements are formed from the rigorous application of energy 
theory to a polynomial representation of displacement functions. An 
important feature in the treatment of bending is that transverse shear 
flexibility is included. 

All of the properties of the elements except those of TRIM6 and TRPLT1 
are assumed uniform over their surfaces. For elements TRIM6 and TRPLT1, the 
thickness can have bilinear variation over their surfaces. In addition, 
element TRIM6 has bilinear variation over the surface for the temperature 
distribution. 

The detailed discussion of the plate elements is divided into subsections 
according to the following topics: membrane triangles, TRMEM, QDMEM, QDMEM1, 

QDMEM2, the basic bending triangle, TRBSC; composite triangles and quadri- 
laterals, TRPLT1, TRIA1, TRIA2, QDPLT, QUAD1, QUAD2, the treatment of inertia 
properties; the isoparametric quadrilateral membrane element, QDMEM1, linear 
strain membrane triangle, TRIM6; higher order bending element, TRPLT1. The 
accuracy of the bending plate elements m various applications is discussed 
m section 15.2, the accuracy of the quadrilateral membrane elements is dis- 
cussed m section 15.3, and the accuracy of the TRIM6 element is discussed in 
section 15.4. 

5.8,6. TRIM6 The Linear Strain Membrane Element 

This element was first formulated by J. H. Argyns and is described m 
references 1 and 2. The present development is based on the derivation m 
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reference 2, and the important characteristics of the element are 
that • 


1. The stresses and strains vary within the element linearly 

2. Bilinear variation m the planar coordinates for the thickness of 
the element is permitted. 

3 Bilinear variation m the planar coordinates for the temperature 
m the element is provided. 

4. Differential stiffness and piecewise linear analysis capability 
are not implemented at present. 

The element is compared for accuracy against theoretical results m 
section 15 4. The calculation of its mass properties is discussed in 
section 5.8.4 

5.8.6. 1 Geometry and Displacement Field 

The geometry of the element is shorn m figure Al. The element has six 
grid points, three at the vertices and three at the midpoints of the sides. 

u and v are components of displacements parallel to the x- and 
y-axes of the local (element) coordinate system The mplane displacements 
at the corners of the element are represented by the vector {u g } where 

{u,} 1 = [ux v x u 2 v 2 u 3 v 3 u 4 v 4 u 5 v 5 u 6 v 6 J (1) 

Let [k ] be the stiffness matrix referred to the vector {u } , 

00 0 

i.e. , 

!XJ ty = {y ( 2 ) 


where the elements of {f g } are the mplane forces at the comers of the 
element. The stiffness matrix [K ge ] is derived by standard finite element 
procedures. 

The u and v displacements are assumed to vary quadratically with 
position on the surface of the element. 
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Figure Al. TRIM6 membrane element m element coordinate system. 
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u = a-i + & 2 X + a. 2 y + a 4 x 2 + a 5 xy + agy 2 (3) 

v = b 7 + bsx + bgy + b^gx 2 + b^xy + b^y 2 (4) 

The quantities aj , a^ , . . . ag , by , bg , . . ^12 niay be 

regarded as generalized coordinates to which the displacements at the corners 

of the element are uniquely related, i.e , the vector of generalized coordinates 
is expressed as 

T 

{a} = [a^aga^agagbybabgbiobnbia J (5) 


In concise form equations (3) and (4) are written as 


6 

a- Z 

1 = 1 


m n 

l i 

a x y 
l J 


12 p* q 

= L b x y 
i=7 


( 6 ) 


m 


For convenience m later calculations, the range of summation is kept as 
1 to 12 for expressions for both u and v , i.e.. 


u 



m n 
i i 

a i x y 


C8) 


V 


12 


i=l 


P i q i 
b i x y 


C9) 


such that 


mi = 0 , 

m 2 = 1 , 

m 3 = 0 , 

m 4 = 2 , 

m 5 = 1 , 

®6 = 0 

ni = 0 , 

n 2 = 0 , 

n 3 = 1 > 

n 4 = 0 , 

H5 = 1 . 

ng = 2 


a i = m^ = m = 0 i = 7 to 12 (10) 
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P7 = 0 , P8 ® l ; P9 = 0 J Pio = 2 ; Pn = l ; P 12 = 0 

q 7 = 0 ; qa = 1 ; q 9 = 0 ; qio = 0 i qil - 1 5 qi2 = 2 

= p^ = q = 0 x = 1 to 6 . (11) 

In matrix notation, the vector {u } is written as 

e 

<V = M {a} C12 > 

where the 12 * 12 £h] matrix can be obtained by substituting the coordinates 
of the six grid points into equations (8) and (9) • Since complete polynomial 
expressions are chosen for the u and v displacements, the inverse of 
M matrix exists. Hence {a} can be expressed as 

{a} = [h] - 1 {u g } (13) 


Bilinear variation in the x- and y-coordinates is assumed for the 
thickness t of the element, i.e., the thickness t of the element at 
any point (x,y) within the element is given by 

t(x,y) = c x + c 2 x + c 3 y (14) 


In concise form, this is written as 


t 


k=l 


k °k 

c ,x y 


(15) 


The thickness of the element at the three vertices is specified as 
ti , t 3 , ts . Hence the coefficients c^ , c 2 , c 3 can be expressed 
as 


t i a + t 3 b 

(a + b) 


( 16 ) 
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t 3 - tl 


C2 a + b 

(17) 

c 3 = “ (t 5 - ci) 

(18) 


where a , b , c are the projected lengths of the triangle on the local 
x- and y-axes and are obtained from the basic coordinates of the vertices 
of the triangle as given m section 4.87.21.2 of the Programmers' Manual. 

The membrane strains are 


£ = = 
x 3x 

12 

£ 

(m -1) n 
l i 

a m.x y 

ii 

(19) 

1=1 



= av _ 
£ y ~ 9y " 

12 

E 

1=1 

P x K' 1 ) 

bVxy 1 
11 J 

(20) 

3u 

Y=: 37 + 

3v 

3x 

/ m x (n -1) (p -1) q v 

£ ( Vi x y + b i p i x y ) 

(21) 


The stress vector {a} is related to the strain vector by the two-dimensional 
elastic modulus matrix, £g ] : 

{a} = [G e ] {e} (22) 

The specification of £g 1 for isotropic and anisotropic materials is the 
same as that given by equations (13), (14), and (15) m section 5.8.4. 

The membrane strain energy of the element is 


E s = j f {ct} T {e} tdxdy 


By virtue of equation (22) and the symmetry of matrix 



(23) 
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( 24 ) 


&.] {e} tdxdy 


Substitution of equation ( 15 ) into equation ( 24 ) results m 


h-kf f» T 

1 1 

<D 

1 1 

t£ > ( 

' 3 r. s. \ 

S v y ) 

' k=l K / 

dxdy 

( 25 ) 

Expressing the elements of 
Gil , G12 , G13 , G 2 2 j 

the symmetric portion of the 
G23 * G33 , i.e.. 

matrix 

C G e ] by 


G ii 

G12 

G 13 






G22 

G 23 



( 26 ) 


sym 


G 33 . 





and performing the matrix multiplication of equation ( 25 ) , the expression for 
strain energy becomes 


E = 


jff { e x 2g U + e / G 22 + Y 2 G 33 + G 12 ( e x £ y + £ y £ x ) 

■I 

V y 7 


+ G13 (b y + YO + G23 ( £ Y + YO 

a a y y 


( 27 ) 



r i V 


dxdy 


To proceed further it is necessary to have a formula for the integral of 
the type 


ff x m y n dxdy 


taken over the area of the element. The value of the integral is given m 
reference 3 as 


14 



m!n! 


(28) 



x™/ 0, dxdy = F(m,n) 

C n-U ( a m*l . ( _ b) m*l 


(m + n + 2) l 


Using equations (19) , 
equation (27) becomes 


(20), (21), and (28) in (27), the first term of 


7 // e x 2Gl1 ( £ c k x ky k ) 


. 12 12 3 

IEEE 

1=1 j=l k=l 


a a c, m m 

i 1 k i j 


dxdy 


F(m x + 


m_ + r,_ - 


(29) 


2 , n + n + s. ) 
i j ir 


Similarly the other terms of equation (27) can be expressed in terms of the 
area integral F . The strain energy, E g , can also be expressed as 


E s -I taf [k gen ](a} (30) 

where E k gen ] the stiffness matrix with respect to generalized coordinates 
{a} Expressing each of the terms of the right-hand side of equation (27) m 
terms of the area integral F and comparing the same with equation (30) , the 
jth element of the ith row of the generalized stiffness matrix is 


k ij = C k [eiim.m.FC^ + + r k - 2, ^ + n^ + s k ) 

+ G 2 2q i q ;) F(p i + + r k , q x + q^ + s R - 2) 

+ G 33 j F (nu + m^ + r fc , ^ + + * k “ 2 ) 

(31) 

* P^/tPi + P, + r ic - 2, q, + + sp 

+ F (m x + Pj + ^ k - 1, n x + q + s fc - 1) 

+ Pi n j F ( m 3 + P x + r k ” n -, + ^i + s k " 1 ^J (continued) 
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* + ?i * r* - 1, n. * q t ♦ s k - 1} 

+ m q.F(m + p. + r, - 1, n + q. + s, -1)1 
i n 3 i k ’ i n j k J 


+ Gi 3 j (nun. + m i n j ) F (n^ + m^ + r k - 1, n ± + + s k - 1) 


C31) 

(concluded) 


s k“ « 


+ m^p i F(m^ + Pj. * r k " 2, n. + q ± + s k ) 

+ nu p^ F (m i + + x k - 2, n ± + q. + s k )j 

+ G 23 | Cp ± qj + P j q i )F(p i + p^ + r k - 1, q x + q^ + 

+ n.q.F(m. + p. + r., n. + q. + s, - 2) 
i n j l k’ i n j k 

+ + Pi + V n j * % + »k - 25 1] 


Using equation (13), the generalized stiffness matrix C^g en l can b® 
transformed to element stiffness matrix £k ^ as 


[kj * C h ''] T l>«J M 


As a final step, the stiffness matrix is transformed from the local 
element coordinate system to the basic coordinate system of the grid points 
and to the global coordinate system. Let the transformation for displacements 
be 


{“basic} ‘ M T {“, 


element 


global 


}' M {“basic} 


K as ic ] ■ M [^element M 
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and 


[“global] ‘ [ T 1 T [“basic] W C363 

Substituting equation (32) m equation (35) and equation (36) , the global 
stiffness matrix becomes 


[“global] - W 1 W [«-‘] T [“ gen ] [H-‘] [*] T [T] (37) 

Equivalent Thermal Load Vector : 

Thermal expansion of an element produces equivalent loads at the grid points. 
Thermal expansion is represented by a vector of thermal strains. 


{e t > = < 

M 

e yt 

> = < 

(a . 
el 

%2 

> (T - T q ) = {« e > (T - T q ) 

(38) 




, “e3 i 




Where {a^ } = [»]-> (a } is a vector of thermal expansion coefficients, 

DO is the strain transformation matrix given m equation (15) m section 
5.8.4, and {ct } is the vector of thermal expansion coefficients m the 
material axis system, T q is the reference or stress-free temperature of the 
material, and T is the temperature at any point (x,y) m the element and 
is given by a linear polynomial. 

T = dj + d 2 x + d 3 y (39) 


In concise form, this is written as 


t„ u„ 


777 V" , 'A & 

T = L d 0 x y 


SL=1 


(40) 
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The temperature Ti , T 3 , and T 5 at the three vertices of the 
element will be modified by the reference temperature T q and used to 
evaluate the three constants di , d .2 , and d 3 : 


* 1 ° • 3 

dl " (a + b) 

(41) 

T’ - T« 
d2 “ (a + b) 

(42) 

d 3 [tj - <g 

(43) 


where 


T i = ( T 1 " V : T 3 = C T 3 - V 5 and T s = (T 3 - T q ) 


(44) 


An equivalent elastic state of stress that will produce the same thermal 
strains is 


<V ■ W {e t } * C G e ] { “e> (¥ - V 


(45) 


An equivalent set of generalized loads {P^^} applied to corners of 
the element is obtained from the relation 


{a)t ^gen 1 = jf {E>t <«t> tdA 


= // (^[Gj { % } (Ej 

( Y' r k s k\ , , 

\ k 1 ° k 7 / 


(46) 
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Performing the matrix multiplications in equation (46) and using the following 
notations, viz 


G H = G H a el + G 12« e2 + Gl 3 a e3 


G 22 = G 12 a el + G 22 a ^9 * G 23 a 


e2 


e3 


G 33 = G 13 a el + G 23 a e2 + G 33 a e 3 


(47) 

(48) 

(49) 


Equation (46) reduces to 


{a}t {P gan } = // tE x G W + E y G i 2 * ^ 


(i vV*) (e 

'£»! 7 v k=l 


r k s k 
c vX y 


dxdy 


(SO) 


Performing the integration term by term, the first term in equation (50) 
becomes 


// Vi, (i vV) (s vV) 


dxdy 


12 3 3 ' 

‘ !?i iS £ ““ViVl x("i * r k * *1 ' 13 

• y^ n i + s k u JlP dxdy 


(51) 


12 3 3 

= 52 52 52 G ii a m.c,d F(m. + r, + t n - 1, n + s, + u„) 

i=l k=I M 11 i i k A 1 k % 1 k 1 


Similarly, the second and third terms of equation (50) reduce to 


12 3 3 

£, £ £ G 22 b i q i C k d » F( Pl ♦ r k + V V - S k * U * ' 15 


1=1 k=i a=i 
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and 


? ? ^ G 3 3 ° kd & { a i n i FCm i 

+ Vx FCp i * r k + h ~ *i + s k + v} 


k + V n i + s k * u £ " ^ 


respectively. From equations (50) and (51) , the ith element of the generalized 

load vector {P } is 
gen 


3 3 


CP 


A * £ £ Vi [g5j» F(m ♦ r fc + t - 1, n x + * k ♦ u > 

k=l £=1 u 

* G 22<li F (Pi * r k * V % * S k * \ - 13 
+ G 33 ( n i Ft \ * r k + V n t * S k + \ - 13 
+ * r k + h ' 1 - \ * S k * U 2 3 } 


(S2) 


The generalized equivalent load vector {P } is transformed to load vector 

gen 

{P e } m element coordinate and to {P } , in global coordinates by the 
following transformations 


{p } = 

e gen 


(53) 


and 


< F g > - M T [e] {P e > 


(54) 


After the grid point displacements have been evaluated, stresses m the 
element are computed by combining the relationships 


<V ■ M M T <y 


(55) 


{a} = £h {u } 


(56) 
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{e} is evaluated from equations (19), (20), and (21). Stress vector 
{o} is then equal to 

{a} = [gJ ({e> - {e t » (S7) 

The stresses are computed at the three vertices and at the centroid. 

The principal stresses and the maximum shear force are computed from the 
elements of {a} . The direction of the maximum principal stress is 
referenced to the side joining grid points 1 and 3 of the triangle. 
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5.8.7 TRPLT1, higher order bending element 

This element was developed by Narayanaswami (refs. 1 and 2) as a 
modification of the high precision bending element of Cowper, et al. (ref. 3). 
The element has grid points at the vertices and at the midpoints of the sides 
of the triangle. At each grid point, there are three degrees of freedom 
in the element coordinate system, viz, the transverse displacement, w , 
normal to the X-Y plane, with positive direction outward from the paper, and 
the rotations of the normal to the plate a and 8 , with positive directions 
following from the right-hand rule. The element, thus, has 18 degrees of 
freedom m the element coordinate system. The transverse displacement, 
w , at any point within and on the boundaries of the element is assumed to 
vary as a quintic polynomial. Since the variation of deflection along any 
edge is a quintic polynomial in the edgewise coordinate, the six coefficients 
of this polynomial are uniquely determined by deflection and edgewise slope 
at the three grid points of the edge. Displacements are thus continuous 
between two elements that have a common edge. The rotation about each edge 
is constrained to vary cubically; however, since the rotations are defined 
only at three points along an edge, there is no rotation continuity between 
two elements that have a common edge. The element thus belongs to the class of 
nonconforming elements. The requirement that the edge rotation varies 
cubically along each edge established three constraint equations between the 
coefficients of the quintic polynomial for w . These equations together 
with the 18 relations between the grid point degrees of freedom and the 
polynomial coefficients serve to evaluate uniquely the 21 coefficients 
aj to a 2 i of the quintic polynomial assumed for v the transverse displacement. 

5. 8.7.1 Derivation of element properties 

Element geometry : Rectangular Cartesian coordinates are used m the 

formulation. An arbitrary triangular element is shown in figure A2. 

X , Y , and Z are the basic coordinates, x , y , and z are the local 
coordinates. The grid points of the element are numbered m counter-clockwise 
direction as shown m the figure. The lengths a , b , and c shown m 
figure A2 can be easily evaluated from the basic coordinate (Xi, Yi, Zi), 

(X3 , y 3 ^ z 3 )> and ( x 5 > y 5 > z s) of the vertices of the triangle. 
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Displacement field : The deflection w(x,y) within the triangular 

element is assumed to vary as a qumtic polynomial in the local coordinates, 
that is. 


w(x,y) = aj + a 2 X + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 + a 7 x 3 
+ a8X 2 y + a 9 xy 2 + aioy 3 + anx 4 + aj 2 x 3 y 
+ a 13 x 2 y 2 + a 14 xy 3 + a 15 y 4 + a 16 x 5 + a 17 x 4 y 
+ a 18 x 3 y 2 + a 19 x 2 y 3 + a 20 xy 4 + a 21 y 5 


In concise form, this is written as 
21 m. n. 

w = * 2 x V 1 

1=1 


(la) 


There are 21 independent coefficients, aj to a 21 . These are evaluated 
by the following procedure. 

The element has 18 degrees of freedom; namely, lateral displacement 

w in the z-direction, rotation a about the x-axis, and rotation 8 

about the y-axis at each of the six grid points. The rotations a and 8 

are obtained from the definitions of transverse shear strains v and 

'xz 

y , that is. 


Y -!*+ B 
’xz 3x 


3 w 

Y = x a 

'yz 3y 


( 2 ) 


It is shown later on that y and y and hence a and 8 at any 

xz yz 7 

grid point can be expressed in terms of the coefficients aj to a 2 i . 

Thus, 18 equations relating w , a and 8 at the grid points to the 21 
coefficients are obtained. Three additional relations are required so that 
the 21 coefficients can be uniquely determined. These relations are obtained 
by imposing the condition that the edge rotation varies cubically along each 
edge. It is clear that these three constraint equations involve only the 
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coefficients of the fifth degree terms in equation Cl) > since the lower 
degree terms satisfy the condition of cubic edge rotation automatically. 
Moreover , the condition depends only on the orientation of an edge. Along 
the edge defined by grid points 1 and 3 (where y = 0) , the condition of 
the cubic edge rotation requires that 

a 17 = 0 (3) 

Along the edge defined by grid points 1 and 5 (me lined at angle 6 to the 

x-axis) the edge rotation r is given by 

0 

r e = 3 sm 6 + a cos 5 = - (5ai 6 x 4 + 4a 17 x 3 y + 3a lg x 2 y 2 

+ 2a! 8 xy 3 + a 20 y 4 ) sin 6 + (a 17 x 4 + 2a 18 x 3 y + 3a 19 x 2 y 2 (4) 

+ 4a 28 xy 3 + Sa 21 y 4 ) cos 6 + . 


where the dots indicate terms of third or lower degree. Also, along this 
edge. 



(S) 


( 6 ) 


By substituting x and y from equation (5) and cos 6 and sm 6 
from equations (6) into equation (4) and rearranging (so that the leading 
terms are positive) , the condition for cubic variation of rotation about 
edge 1-5 is 


5b 4 ca 16 + (4b 3 c 2 - b s )a 17 + (3b 2 c 3 - 2b 4 c)a 18 

+ (2bc 4 - 3b 3 c 2 )ai9 + (c s - 4b 2 c 3 )a 28 - 5bc 4 a 21 = 0 
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Similarly, the condition for cubic variation of the rotation about the edge 
defined by grid points 3 and 5 (fig. 1) can be written as 

Sa^ca^g + (-4a 3 c 2 + a 5 ) a 17 + (3a 2 c 3 - 2a if c)a 1 g 

( 8 ) 

+ (-2ac lf + 3a 3 c 2 )a lg + (c 5 - 4a 2 c 3 )a 20 + 5ac**a 2 i = 0 


The 18 relations between grid point displacements and the coefficients of the 
polynomial in equation (6) are written as 


{5} = [Q] {a} 


(9) 


where {6} is the vector of grid point displacements, [Q] is the 

18 x 21 matrix involving the coordinates of grid points substituted into the 

functions w £eq. (1)J and the appropriate expressions of a and g derived 

m detail later, and {a} is the column vector of coefficients a j to 

a 2 i * The [q] matrix is now augmented by the three constraint equations 

(3) , (7) , and (8) to form a new 21 * 21 matrix M in the following equation* 

{6 a > = [R] {a} (10) 


where 


« a >- < 


(U}\ 
o 
0 

x 0 , 


(10a) 


For use in the evaluation 
m terms of {$ } ; and, 

cl 

matrix exists. The 

T-15 and T-21 elements 
the polynomials for w . 
(ref. 3) give an explicit 


% 

of the stiffness matrix, {a} needs to be expressed 
hence, it has to be established that the inverse of 
nonsingulanty of such a matrix M for the 
of Bell (ref. 4) follows from the completeness of 
For the high precision element, Cowper et al. 
expression for the determinant of such a matrix and 
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show that the matrix is nonsingular m all practical situations. For this 
element, a numerical experiment described m reference 1 verifies that 
R is nonsingular for all practical cases. Hence, equation (10) is inverted 
to give 


{a} = [r]" 1 {6 } 

cl 


(ID 


This equation can also be written as 

{a} = [s] {6} (12) 

where []s]j is a 21 x 18 matrix and consists of the first 18 columns of 

DO - 1 • 

From the computational standpoint, it is advantageous to substitute 
equation (3) into equation (1) and replace coefficients a lg to a 21 
by coefficients a 17 to a 20 > respectively. The matrix £q[] then is of 
size 18 x 20, £r] becomes 20 x 20 , and [s] becomes 20 x 18. To add 

to the clarity of presentation, however, the complete qumtic polynomial for 
w m equation (1) is retained throughout this section, and matrices £q[] , 
[R] , and [s] and vector {a} will have sizes 18 x 21, 21 x 21, 21 x 18, 
and 21 x l, respectively. 

Elastic relationships The elastic relationships are obtained from the 
theory of deformation for plates (ref 5) The curvatures are defined by 


/ X ^ 


_ 3|!_ 

\ 

X 


3x 


X 

> < 

3a 

> 

y 

3y 




3a 

M 

xy J 


3x " 

3y J 


Bending and twisting moments are related to curvatures by 


(13) 
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( 



= [D] 


xy t 


(14) 


where M is, in general, a full symmetric matrix of elastic coefficients. 
For a solid isotropic plate of uniform thickness t , 


M 


Et 3 

12(1 - v 2 ) 



(15) 


The thickness t of the element is assumed to vary bilinearly with 
position over the surface 

t = ci + C 2 X + C3y (16) 


In concise form, it is written as 


t = 



C k X 


T k. S k 


(16a) 


The thickness of the three vertices of the element ti , t 3 
will be used to evaluate the constants Ci , C 2 > and c 3 
shown that 


and t 5 
It can be 


tia + t 3 b 
Cl (a + b) 

(17) 

*3 - ti 

C2 “ (a + b) 

CIS) 

c 3 = ^ [ts " c lJ 

(19) 
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where a , b , c are the length of the element marked m figure Bl. 
For an isotropic plate, [d] becomes 


M ■ h ^ ( £ £ £ 

x 1=1 1=1 k=l 


r i* r , +r k s i +s / s k 

c i c : c k x y 


where 


1 



[G 


e 


] = 




0 

0 


0 


0 


E 

2(1 + v) J 


( 20 ) 


( 21 ) 


For anisotropic materials vrith the material orientation axis inclined at <f> 

to the x-axis, the material elastic modulus matrix [d 1 is transformed to the 

L m J 

element elastic modulus matrix by 

[D] = [U] T [D m ] [U] (22) 


inhere 


[U] = 


cos 2 <j> 
sm 2 <j) 

-2cos <p sm <j> 


sm 2 <|i 
cos 2 $ 

2cos 4> sin ^ 


cos <j> sm <ji 
-cos (j> sin <Ji 
cos 2 (j> - sm 2 $ 


(23) 


Tile positive sense of bending and twisting moments and transverse shear 
resultants is shown m figure A3. 

The moment equilibrium equations are written as 


V + 
x 



3M 


xy 


3y 


= 0 


(24) 
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0 


(25) 


V + 

y 


3M 3M 
y . xy _ 

3y 3x 


Transverse shear strains are related to the shear resultants by 






(26) 


The matrix [j] is, m general, a full symmetric 2x2 matrix of elements 
^11 > J 12 (J 21 = J 12 ) and J 22 • P° r a plate with isotropic transverse 
shear material. 


M = 


1 

Gt* 


1 0 
0 1 . 


(27) 


where G is the shear modulus and t* is an "effective" thickness for 
transverse shear. For a simple case of a plate of uniform thickness t , 
t* has the value t . 

From equations (24), (25), and (26), it follows that 


(28) 


Performing the partial differentiation with respect to x and y on 
equation (14) , with subscripts on D denoting the elements of [d] , 
results m 




3M 

3M 1 


r3M 

3M 1 

Y xz = 

-Jn 


-J3L 

3y _ 

- J12 

U- 

xy 

3x 



raw 

3M ' 


r 3M 

3M 1 

Y yz = 

-J12 


xy 

3 y . 

- J22 

i 

a? 

05 

1 
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3M 

x _ 

3x " 



3y 

Di2~g^ + D13 


3x xy 

3x 


\ 


X 3x x 9x y 9x xv 

17 ” ° 12 17 + D22 1^ + D 23-a^ 


3M 


3X x 9 X^ 

Dl3 ~8x + D23 “fl^ * ° 33 ' 


9 X 


2QL 


3x 


3x 


3x 


3M xy 3x x 3x y 3x 

= °13^y + D 2 3-g^ + D 33 -^ 


xy 


( 29 ) 


where the symmetry of the M matrix has been used. By substituting 
equations (29) into equations (28) , 


Y = 
xz 


= - Ji 


F 3x x 3x y 3x xy 

1 L Dll_ ^ + Dl2 ~l7 * Dl3_ 9x 

3X x 3X„ 8X^1 

+ d i 3- 97 + D23 i7 + D33 "§t : J 

[■ 




9 X 


J 1Z d 12‘ 


3 X 


ly + ° 22 ~dp + ° 23 3y 


221 


3X, 


3 X V 


+ Dl3 ~57 + ° 23 ^ + ° 33 "iF 


(30) 


and 




9 X. 


9 x. 


+ D 13^7 + D23 ~^ 


r 3x x 

- J22[D 12 - 1 y 

9 X. 


>x y 


° 22 — 


+ ° 13 ^7 + ° 23 ‘ 


( 31 ) 
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Rearranging and writing equations (30) and (31) m matrix notation yields 

/ \ 


M . 

A 11 A 12 A 13 A 14 A 15 A 16 

U ' 

_ A 2 1 A 2 2 A 2 3 A 24 A 2 5 A 2 6 


X x,x 

X y,x 


I X xy,x I 
I X x,y I 


x 


y,r 


x 


xy,y 


(32) 


\ / 

where a comma m the subscript denotes partial differentiation and where 


An = “( J 11 D 11 * J 12°13) (33a) 
A 12 = “( J 11 d 12 + J 12 D 23) (33b) 
A 13 = -C J llDi3 + Jl2 D 33) (33 c) 
a 14 = “(JllDl3 + J 12°12) (33d) 
A 15 = _ ( J 11 d 23 + j 12 D 22) (33e) 
A 1 6 = “( J 11 D 33 + J 12°23) (33f) 
a 21 “ -(Jl2°ll + J22°13) (33g) 
A 22 ~ “( j 12 D 13 + J 22 D 233 (33h) 
A 23 = -(Jl2 D 13 + ^22°33) (33l) 
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A 2 4 = -(J12D13 + ^22^12) 


(33 j) 


A25 = -(J12D23 + J22°22) 


A 26 ~ -(J 12 D 33 + J 22°23) 


From equations (2) and (13), it follows that 


38 9£w 3y xz 

Xx = ’ 3x = 9x 2 " 3X 

3cc 9 2 w 3y yz 

V " 37 = 3y 2 ” 37 



9a 36 „ 2 d 2 w 3y xz 
9x " 9y - ^ 9x3y ~ 9y 



/ 


Shear forces (and hence shear strains) are proportional to the third 
derivatives of the displacements. Since the displacement within the 
element is assumed to vary as a quintic polynomial, shear strains are 
expressed by a quadratic polynomial as follows 

Y^„ = t»i + b 2 X + b 3 y + bi^x 2 + b 5 xy + b 6 y 2 


Y v _ = b 7 + b 8 x + b 9 y + bi 0 x 2 + b lx xy + b 12 y 2 

The task now is to express the unknown coefficients b x to b 6 and b 
b 12 m terms of the generalized coordinates a x to a 2X . Performing 
differentiations on x x > Xy > and and substituting w , y^ , 

and y from equations (1) , (35) , and (36) into equations (34) 


(33k) 

(331) 


(34) 


(35) 

(36) 

- to 
the 


34 



= 6a7 + 24a 11 x + 6a 12 y + 60a 16 x 2 


3 d w 


3 2 Y 


xz 


X > X 3x 3 3x 2 

+ 24ai7xy + 6ai 8 y 2 - 2bi 


( 37 ) 


3 3 w 


3 2 Y. 


3x3y* 

+ 12a ig xy + 12a 20 y 2 - b u 


— = 2ag + 4ax3X + 6aj 4 y + 6a 18 x 2 


(38) 


V.* 2 a „2 a „ 3x3y 


3\ *\z !!v = 4aa * 12ai2X * 8aisy 


3x 2 3y 3x 

+ 24a^x 2 + 24a^gxy + 12a ig y 2 - b 5 - 2bio 


(39) 


3 3 w 


3 2 Y 


XZ 


X ’ y 3x 2 3y 3x3y 

+ 12a 18 xy + 6a 19 y 2 - b 5 


= 2a 8 + 6a^ 2 x + 4a ig y + 12a! 7X 2 


(40) 


3 3 w 


3 2 Y 


yz _ 


™ 9y 3 9y2 

+ 24a 28 xy + 60a 2 iy 2 - 2b^ 2 


= 6a^o + 6a 14 x + 24a 15 y + 6a 19 x 2 


( 41 ) 


and 


X - 2 
xy,y 


3 z w 


3 2 Y 


xz 


3 2 y 


y z _ 


3x3y 


3x3y 2 3y 2 

+ 12ai 8 x 2 + 24ai 9 xy + 24a 20 y 


4ag + 8ax3X + 12a xi*y 
2 - 2b 6 - b n 


(42) 


By substituting equations (35) to (42) into equations (32) , the following 
equations are obtained. 
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bi + b2X-+ b3y + b^x 2 + bsxy + b6y 2 = An(6a7 + 24anX 
+ 6ai2y + 60aigx 2 + 24anxy + 6ai 8 y 2 - 2bi f ) 

+ Ai 2 (2a9 + 4ai 3 x + 6aii*y + 6ai 8 x 2 + 12a 19 xy + l2a 2 oy 2 - bn) 

+ A 13 (4a 8 + 12a 12 x + 8aj 3 y + 24a! 7X 2 + 24ai 8 xy + 12aigy 2 

(43) 

- b 5 - 2b 10 ) + A lt} .(2A 8 + 6a 12 x + 4a 13 y + 12a 17 x 2 + 12a 18 xy 

_ \ , 

+ 6aigy z - bs) + Ai5,(6aio + 6ai4X + 24ajsx + 6ax 9 x^ 

+ 24a 28 xy + 60a2iy 2 - 2b^ 2 ) + Ai 8 (4ag + 8a 13 x + 12ai4y 
+ 12ai 8 x 2 + 24aigxy + 24a 28 y 2 - 2b 8 - bn) 

b7 + b 8 x + bgy + bi 0 x 2 + bnxy + b^x 2 = A 21 (6a7 + 24anx 
+ 6ai 2 y + 60a x 8 x 2 + 24a*7xy + 6ai 8 y 2 - 2b 4) 

+ A 22 (2ag + 4a i3 x + ba^y + 6aj 8 x 2 + 12aigxy + 12a 20 y 2 

- bn) + A 23 (4a 8 + 12a 12 x + 8a 13 y + 24a i7 x 2 + 24a i8 xy 

(44) 

+ 12aigy 2 - b 5 - 2bio) + A 24 (2a 8 + 6a 12 x + 4a 13 y + 12a 17 x z 
+ 12aj 8 xy + 6a 19 y 2 - bs) + A 2 s(6aio + 6a X 4X + 24aisy 
+ 6aigx 2 + 24a 28 xy + 60a 2 iy 2 - 2b i 2 ) + A 28 (4a 9 + 8ai 3 x 

+ 12ax4y + 12ai 8 x 2 + 24ajgxy + 24a 28 y 2 - 2b 8 - bn) 

By comparing coefficients of like powers m x , y , x 2 , xy , and 
y 2 and constants of equations (43) and (44), the coefficients b^ to 
b 6 and b7 to b 12 can be expressed in terms of the generalized 
coordinates aj to a 2 j . Thus 
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b2 = 24Anan + 6(Am + 2Ai3)ax2 + 4(Ai2 + 2Axg)ax3 + 6Axsax4 
b3 = 6Anai2 + 4(Ai4 + 2Ai3)ax3 + 6(Aj2 + 2Axg)ax4 + 24Axsaxs 
b4 = 60Aiiai$ + 12(Ai4 + 2Ai3)ai7 + 6CAx 2 + 2Aig)a-18 + ^Axsaxg 
bs = 24Axx a x7 + l2(Aj4 + 2Ax3)axs + 12(Aj2 + 2Axg)axg + 24Axsa 2 Q „ 
bg = 6Anaig + 6(Ai4 + 2Ai3)a ig + 12(A 12 + 2Ajg)a 2 o + 60Axsa 2 x 


> (45) 


bi = 6Ana7 + 2(Ai4 + 2Ax3) a 8 + 2 (Aj 2 + 2Aig)ag + 6Ax5axo 

- 2Ajib4 - (Ax 3 + Ai4)b5 - 2 A^gbg - 2 Ax 3 bxo - (Ax2 + Ax 6 )bll 

- 2A]i5b 12 


bs = 24A2x a n + 6(A 2 4 + 2A 2 3)aj 2 + 4(A 22 + 2A 2 g)ax3 + 6A 25 ax4 
bg = 6A 2 x a x 2 + 4(A 2 4 + 2A 2 3)ax3 + 6(A 22 + 2A 2 g)ax4 + 24A 25 ai 5 
bjo = 60A 2 iai6 + 12(A 2 4 + 2A 2 3)ax7 + 6(A 2 2 + 2A 2 g)ai8 + 6A 2 saxg 


bn = 24A 2 jax7 + 12 (A 2 l + 2A 2 3)axg + 12(A 22 + 2A 2 g)axg + 24A 2 sa 2 o 


> (46) 


bi 2 = 6A 2 iai8 + 6(A 2 4 + 2A 2 3)a ig + 12(A 22 + 2A 2 g)a 2 g + 60A 2 sa 21 


b 7 = 6A 2 x a 7 + 2(A 2 4 + 2A 23 )a8 + 2(A 22 + 2A 2 g)a g + 6A 2 gaxo 

- 2A 2 ib4 = (A 2 3 + A 2 4)bs - 2A 2 gbg - 2A 2 3bxo - ( A 22 + A 2 g)bxx 

- 2A 2 5bi 2 
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If equations (45) and (46) are substituted into equations (35) and (36) , 
the explicit relation between the transverse shear strain and the generalized 
coordinates (i. e. , coefficients of the displacement polynomial) can be 
obtained m matrix notation as 

iy} - [Bi] {a} (47) 

where [Bi] is a 2 x 21 matrix whose nonzero elements are as follows: 


B l(l,7) = 6 Au 

(47a) 

Bid, 8 ) = 2A 31 

(47b) 

Bi(l,9) = 2A 32 

(47c) 

Bx (1,10) - 6Axs 

(47d) 

Bx (1,11) = 24A xl x 

(47e) 

Bx(l»12) = 6 (A 3i x + Axxy) 

(47f) 

Bi(l,13) = 4 (A 32 x + A 31 y) 

(47g) 

Bx(l,14) = 6 (Ax 5 X + A 32 y) 

(47h) 

Bx (1,15) = 24A 15 y 

(47i) 

Bx (1,16) — -120(A 2 x x + Ax 3 A 2 x - 0.5AxxX 2 ) 

(47 j) 


38 


Bi(1 3 17) - -24 [Aj 1 (A 31 + A 38 ) + A 13 A 33 + A 21 A 3 g 
- 0.5A 31 x 2 - A n xy] 


(47k) 



= - 12 (A 11 A 32 + A13A34 + ^38^31 + A39A33 + Aj^Ajg 

C 471 ) 

+ A15A21 - 0 . 5 A 32 x 2 - A 31 xy - 0 SA n y 2 ) 

61(1,19) = -12(A 11 A 15 + A 13 A 2 5 + A 3 qA 3 2 + A 3 gA 3 4 + A 16 A 31 

(47m) 

+ A15A33 - 0 5 A 15 x 2 - A 32 xy - O.SA 31 y 2 ) 

B l(lj 20 ) = - 24 CA 15 A 38 + A 25 A 39 + AigA 3 2 + A 15 A 34 


- Ai 5 xy - 0.5A 32 y 2 ) 

(47n) 

B 1 ( 1 , 21 ) = - 120 CAi 5 Aig + A 15 A 25 - 0 5Ai 5 y 2 ) 

(47o) 

63 .( 2,73 = 6 A 21 

(47p) 

b i(2,8) = 2A33 

(47q) 

(2,9) = 2 A 3 4 

(47r) 

B 1 (2,10) = 6 A 25 

(47s) 

Bi (2 ,11) = 24A 2i x 

(47t) 

b i (2,12) = 6 (A 33 x + A 2 i y) 

(47u) 

B i (2,13) = 4 (A 34 X + A 33 y) 

(47v) 

Bi (2,14) = 6 (A 25 x + A 3 £f y) 

(47w) 

Bi (2,15) = 24A 25 y 

(47x) 
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Bi(2,16) = -120 (AijA21 + A 23 A 21 - O. 5 A 21 X 2 ) 


C47y) 


where 
A 25 > 


BiC2,17) = -24(A 2 iA 31 + AnAtj .0 + A 23 A 33 + A 21 A 34 

- 0.5A 33 x 2 - A 2 ].xy) 

Bi (2,18) = -I 2 CA 21 A 32 + A 23 A 34 + A 40 A 31 + A 41 A 33 + A 26 A 11 
+ A 25 A 21 - O. 5 A 34 X 2 - A 33 xy - 0.5A 21 y 2 ) 

Bi (2,19) = -12(A 2 iAi5 + A 23 A 25 + A 40 A 32 + A 42 A 34 + A26A 31 
+ A 25 A 33 - O. 5 A 25 X 2 - A 3 4 xy - 0.5A 33 y 2 ) 

B l(2>20) = -24(Aj5A4o + A 25 A 41 + A 26 A 32 + A 25 A 34 - A 2 sxy 

- 0. 5A 34 y 2 ) 

B 1 (2,21) = -120CA 15 A 26 + A 2 25 - O.SA 25 y 2 ) 

A 11 > A 12 > A 13 > A 14 > A 15 > Aig , A 2 l , A 22 > A 23 
and A 2 g are as defined in equations (33) and 

\ 

A 31 = A X 4 + 2A 13 

a 32 = A 12 + 2A 16 

A 33 = A 2 4 + 2 A 2 3 V 

A34 = A 2 2 + 2A 2 g 

A 35 " a 33 + A 11 


(47z) 

(47aa) 

(47bb) 

(47cc) 

(47dd) 

, A 24 

(48) 

(continued) 
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A 36 “ A 34 + A 31 


a 37 = A 25 + a 32 
A 38 = A 1 3 + A 14 
A 39 = A 12 + A 16 
a 40 = A 23 + A 24 
A 41 = A 22 + A 26 


(48) 

(concluded) 


If the plate is assumed to be rigid m transverse shear, the coefficients 
An to a 16 and A 2 l to A26 of equations (33) are zero (since G = «) , 
and, hence, coefficients bj > to b 6 and b 7 to bi 2 of equations (40) 
and (41) are zero. Moreover, the transverse shear strains vary linearly with 
G — 1 with {y> approaching 0 as G -*• » , that is, convergence to the 
limiting case of zero transverse shear is uniform. 

Stiffness matrix . The strain energy for a plate may be written as 

u = |//(( m}T (X> + (V} T (Y>) dxdy (49) 

where {M} is the vector of bending and twisting moments per unit length, 

{y} is the vector of curvatures, {V} is the vector of transverse shear 
forces per unit length, and {y} is the vector of transverse shear strains. 
Substituting equations (14) and (26) into equation (49) , and using the 
symmetry of [d] and [j] matrices, yields 

u = jffW 7 [D] *X> + (T> T [G] (y> dxdy (SO) 


where 
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[G] = [J]- 1 


With en l denoting the generalized stiffness matrix, that is the 
stiffness matrix with respect to generalized coordinates (coefficients of 
the displacement polynomial) {a} , the strain energy can also be expressed 


u - i <a> T [K gen ] {a} 


The vector of curvatures {y} is now rewritten as 
- <X1> * tx 2 > - ([Ba] + [B S ]) {a} 


(m. -2) n. 

Zj a i m i (nu - l)x y 


m (n -2) 

Ea.n^n. - l)xV 1 


(m.-l) (n -1) 

£ a i m i n i x y 


{X2> = 


(53b) 
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It follows that (xi) is the vector of curvature m the absence of 
transverse shear and (x 2 } th e contribution of transverse shear to the 
vector of curvatures. 

Substituting equation (53) into (50) and comparing the resultant equation 
with (52) , noting that {a} is independent of x and y , the generalized 
stiffness matrix can be obtained as 


C^g en 3 = ff [B 2 f M [B 2 ] dxdy +// [b 2 ] T [d] [b 3 ] dxdy 

* ff 1>3] T D>] [B 2 ] dxdy + ff [B 3 ] T [D] [b 3 ] dxdy (54) 
+ ff [Bl] T [G] [Bj dxdy 


The evaluation of the elements of the generalized stiffness matrix 
[k 1 m closed form is, though straightforward, very tedious The 
first term J j [B 2 J [d] [PlA dxdy is evaluated m closed form, the other 
four terms are evaluated by using numerical integration. If the plate is 
assumed to be rigid m transverse shear, the matrices [B^] and [b 3 ] are 
null, and the last four terms vanish. The numerical integration formulae 
used are the seven-point integration scheme listed m reference 6 and are 
given below for easy reference. For a triangle, the integrals of the form 


1 1-Li 

/ f 

0 0 




(55) 


can be integrated by using a seven-point numerical integration which can 
exactly integrate functions up to and including qumtic order The value 
of the integral is given by 


i= E 

k=l 


W. 




(56) 
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where the points and the weighting factors are as follows 


Point 

Triangular Coordinates 

Weight, 2W^ 


hi. 


l 3 


1 

1/3, 

1/3, 

1/3 

0.225 

2 

a l 

3i 

3i 


3 

t 

4 

3l 

“1 

3i 

0.13239415 

3l 


“1 


5 

<*2 

32 



6 

32 

a 2 

32 

0.12593918 

7 

32 

32 

a 2 




with 


ct! = 0.05971588 8j = 0.47014206 

ct 2 = 0.79742699 8 2 = 0.101286505 

Note the error m the value of as given m reference 6, page 151. 

Denoting by Gjj , Gj 2 , Gj3 , G 22 , G 2 3 , and G33 the symmetric 

portion of the G matrix of equation C21) , it can be shown that the jth 
element of the ith row of the generalized stiffness matrix IXeJ , for the 
case of a plate infinitely rigid m transverse shear, is given by 
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(K J 


1 J ' gen 


, .5 s o 

= X y y y c . c, c 

12 k * k * 


[G u m m ^ - 1)0^ - D F^ + ^ + r ki - * k . 


+ r, - 4, n + n + s, + s, + s, ) 
k 3 i J k l k 2 k 3 


+ - 1)0^ - 1) + m. + 

+ r, + r. , n + n + s v + s, + s. - 4) 

k 2 k 3 i 3 k l k 2 k 3 

+ ( 4G 33 m i m j n i n j + (n^ ~ *) ( n -, ~ *) 

+ mn^ - l)Cn x - 1)}) ® 

+ r, + r. - 2. n + n + s, + s. + s, - 2 ) 

k 2 k 3 1 J k l k 2 k 3 

+ 2G 13 {m i m^n^ (m^ - 1) + nun^ (m^ - 1]} F(>\ + 


+ r kx + r k 2 + r k 3 “ 3 ' n i + n D + S ki 


+ S k 2 + S k 3 ’ 1} + 2G2 3 im J n i 11 ] Cn i ' 1} 

+ m A n ] (n ] ‘ 1} i FCm i + m j * % + r k 2 

+ r k 3 - X > n i + n j + % + % + s k 3 ~ 3) ] 


(57) 


All computations involved m evaluating C^ oen l ^ or tke case a 
plate infinitely rigid m transverse shear can be carried out within the 
computer. For plates with transverse shear flexibility, the contribution 
of the last four integrals of equation (54) will be evaluated using the 
numerical integration formula [eq. (56)] and algebraically added on to the 
closed form expression for C k g en 3 evaluated by equation (57) . 
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Once the generalized stiffness matrix is evaluated, the 

element stiffness matrix in the local element coordinates [K ee 3 is 
obtained as, by virtue of equation (12), 

[K ee ] = [S] T [K gl J [S] (58) 

CK ge l can then be transformed to the global coordinate system of the 
surrounding grid points m the same manner as for all other elements. 

Let the transformation for displacements be 

“basic * M T “element ^ 


and 


“global * M “basic (60) 

Then, the stiffness matrix in global coordinates is 

Mglobal - W T W &] T 0] (61) 

Equivalent thermal bending loads 

The following derivation to obtain the equivalent thermal bending loads 
is given for the case of different thermal gradients at the three vertices 
of the element. This capability is not currently operational m NASTRAN. 
However, the derivation is valid for cases with the same thermal gradients 
at the vertices, if T3 and T5 in equations (83), (74), and (75) are set 
equal to Tj . 

The stress-free strains developed m a free plate due to a variation 
of temperature with depth are 
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(62) 



< T - T raf> * ( “a> < T ' W 


where T is the temperature at any point (x,y,z) of the element, ^ re £ 1S 

the reference or stress-free temperature of the material, and {a } is the 

0 

vector of thermal expansion coefficients m the element coordinate system 
An applied stress vector which would produce the thermal strains is 


V = ^ {e t> = <“a> < T - T „fP 


(63) 


where eg is the matrix of elastic coefficients at the point on the 
cross section 

The generalized equivalent thermal load vector ^Pg en ^ 1S obtained 
as 


{P gan } = >51 / {e>T {< ’t Wv 

The strains (e) are related to the curvatures {y} by 

{e} = -z{y> (65) 

where z is measured from the neutral surface of the plate Substituting 
equations (63) and (65) into equation (64) , 

{P gen } - - W / *<X> T [«.]<«.> ’ W dV 

The variation over the surface of the element of the mean temperature, 

T q , and the thermal gradient at a cross section, T* , is assumed to be 
a bilinear polynomial. 
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T = 
o 



C67) 


T* 



A 1 

d.x y 


( 68 ) 


so that the temperature at any point (x,y,z) T , is 


T = T + T'z (69) 

o 

The constants d. and d.' are evaluated from the values at the 

i l 

vertices . Thus , 

Toia + Tj 3 b 

d l = (a + b) ■ < 70 > 

T03 “ Toi 

dz = (a + b) (71 ^ 

d 3 " ^ [ T 05 - (72) 




( 75 ) 


where Tgj , Tq 3 , and Tq 5 are the difference between the grid point 
mean temperature T 0 \ > T 03 and T os , at grid points 1, 3, and 5, 
respectively, and the reference temperature, T re ^ . 
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It is convenient to define the equivalent thermal moment vector 


<»y - - / [<g v o' - \ ef ) 


zdz 


+t/2 

/ [G J (a } (T’ + T 1 ) zdz 

-t/2 e e o 


( 76 ) 


[<g (« e > t- 


Substituting for t from equation (16a) and for T 1 from equation (68), 


{M t } = - 


12 {c e } T ' 2 E E E c 11 c l2 c l3 

ii=l i 2 =l 13=1 J =1 3 


(r +r +r +p ) (s +s +s +q ) 

il 12 13 V il 12 13 r 


(77) 


At the three vertices the value of {M^} will be given by 


{M t h 

{ M t ) 3 

(M t } 5 


tl 

[G e ] (« e > iJ T{ 

(f 


r 1 ^ 

<“e } 12 T 3 


ts 

^ <“e } 12 T S 


(78) 


(79) 


(80) 


where ti , t 3 , and ts are the thicknesses at the vertices Gi , 

G 3 , and G$ , respectively, of the element. 

I 

The "effective thermal gradient," T , at the vertices is defined as 


T i " i / T i 1 dz 


(78a) 
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Tl * Ti / Ts 1 dz 

(79a) 

t 5 - J x f T 5 z dz 

(80a) 


This capability of specifying the thermal gradients or the thermal 
moments at the three vertices of the element is not currently implemented. 

The theoretical derivations of the evaluation of the thermal load vector 
is, however, given for such linear variation of the thermal gradient values 
over the surface of the element. 

Substituting equations (16), (53), (67), (68), and (69) into equation (66). 

(P gan } ' ' EfW //(<«> * <X*>) T f“ e } 


3 3 3 3 


E E E E 

il=l i 2 =l i3=l 3=1 


c 


c. c. 
*2 x 3 


d« 

3 


(81) 


x 


(r +r. +r +p ) (s. +s. +s. +q.) 

il i 2 13 3 11 12 13 3 

y dxdy 


As in the case of the derivation of the generalized stiffness matrix, the 
generalized thermal load vector will be evaluated m two stages, viz., the 
closed form expression [ P g en li due to [xi] 5 the vector of curvatures 
in the absence of transverse shear, and the numerically integrated 
expression [P^ en l 2 due to [y 2 ]| , the contribution of transverse 

sheat to the vector of curvatures. Using the following notations, viz.. 


G ll " Gll<2t e 1 + Gl2a e 2 + Gl3<x e 3 ( - 82 - ) 

G 22 = G 12 a ei + G 22«g 2 * G 230 es (83) 

G 33 = Gl3 ° t e 1 + G2 3« e2 + G 33« e3 (84) 
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the ith element of the generalized load vector {P x } will be given by 

& 

(P een !l " TJ £ £ £ £ c 1 , c 12 c 1 , d : 

8 1 11=1 i 2 =l l 3 =l 1=1 1 2 3 ] 

' « F(m i * r ij + r i 2 + r i 3 + P J ‘ 2 ' 

"l * s i, * s l 2 + s i 3 * 'll 5 * G 22 n i< n i - D 

Ffm+r + r +r + p , n + s + s 
1 11 1 2 13 *J 1 11 12 

+ s + q - 2) + GI 0 m.n F(m + r + r. 

13 n J 1 33 1 1 v 1 i x i 2 

+ r +p - 1, n. + s + s. + s + q - 1)1 

13 3 l 11 12 13 3 J 


(85) 


The load vector { <ren^ 2 1S ©valu^ed using numerical integration 

and fp^ 1 is obtained as the sum of Ep?Lt.1i and [p^ 1? For plates 

gen J § en 1 gen JZ r 

infinitely rigid m transverse shear, [p^ ] 2 is null. The equivalent 

gen 

thermal bending load m the local element coordinate system is obtained 
as, by virtue of equation (17), 


{P^} = [S] T {P* } 

e u gen 


( 86 ) 


The load vector can then be transformed to the global system by 


(P b = (T] T [E] {P*} (87) 

to 6 

/ 

RECOVERY OF INTERNAL FORCES 

The bending moments and shear forces are recovered at the three 

vertices, the stresses are evaluated at the three vertices and at the centroid 

of the element. After the displacements of the element are transformed 

from the global system {u ] to the element coordinate system £u , 

8 © 
the generalized coordinates {a} are evaluated from equation (12) The 

curvatures (y) are evaluated from equation (53) with the nonzero elements 

of [B 3 ] being as listed below. 
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B 3 C1,11) = -24A n . 


83(1,12) = -6A 31 
63(1,13) = -4A 32 
B 3 (1,14) = -6A 15 
B 3 (1,16) = -120A n x 
B 3 (1,17) = -24(A 31 x + A n y) 
B 3 (l,18) = -12(A 32 x + A 31 y) 
B 3 (l,19) = -12(Ai 5 x + A 32 y) 
B 3 (l,20) = -24A i5 y 
B 3 (2,12) = -6A 21 
B 3 (2,13) = -4A 33 
B 3 (2,14) = -6A 34 
B 3 (2,15) = -24A 25 
B 3 (2, 17) = -24A 21 x 
B 3 (2,18) = -12(A 33 x + A 21 y) 
B 3 (2, 19) = -12(A 34 x + A 33 y) 
B 3 (2,20) = -24(A 25 x + A 34 y) 
B 3 (2,21) = -120A 25 y 
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B 3 (3,ll) = -24A 21 

B 3 (3,12) = - 6 (A U + A 33 ) 

B 3 (3,13) = -4(A 31 + A 34 ) 

B 3 (3,14) = - 6 (A 32 + A 25 ) 

63(3,15) = -120A 15 

^ B 3 (3,16) as -120A 21 x 

B 3 (3, 17) = -24 [(A n + A 33 )x + A 2 iy] 

B 3 (3,18) = -12 [(A 34 + A 31 )x + (A34 + A u )y] 

B 3 (3,19) = -12 [(A 2 5 + A 32 )x + (A 3 4 + A 31 )y] 

B 3 (3,20) = -24A 15 x + (A 32 + A 25 )y 
B 3 (3,21) = -120A 15 y 

inhere An , Ai 2 , , A 34 are as given m equations (33) and (48) 

Moments at the vertices are then obtained from 

(M>i = [d] x { X > - {M t h ( 88 ) 

{M } 3 = [d ] 3 { X > - {M t ) 3 (89) 

Ms = Ld ] 5 - {M t ) 5 (90) 
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The transverse shears are evaluated as follows: 

{y} is evaluated from equations (28) and (47). 

{V} is then evaluated from equations (24) and (25) . 

The stresses at the three vertices are evaluated at distances Zll, 

Z21, Z13, Z23, Z15, and Z25 specified by the user. The stresses at the 

centroid are evaluated at the top and bottom fibers of the element. 
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5.14 TRSHL, Higher Order Shallow Shell Element 

This element was developed by Narayanaswami (ref. 1) . The element has 
grid points at the vertices and at the midpoints of the sides of the triangle 
At each grid point, there are five degrees of freedom m the element coordinate 
system* viz., the membrane displacements u and v parallel to the x and 
y axes, the transverse displacement, w , m the z-direction normal to the 
x-y plane, with positive direction outward from the paper, and the rotations 
of the normal to the shell a and 3 , about the xz and yz planes, with 
positive directions following from the right-hand rule. The element, thus, has 
30 degrees of freedom m the element coordinate system. 

The membrane displacements u and v for the shell are expressed as 
quadratic polynomials and are the same as for the higher order membrane 
triangular element, TRIM6. The displacement function for the normal deflection, 
w , is taken as a qumtic polynomial as for the higher order bending triangular 
element, TRPLT1 . The geometry of the shell surface is approximated by a 
quadratic polynomial m base coordinates Shallow shell theory of Novozhilov 
(ref. 2) is used for including the membrane bending coupling effects. Thus, the 
element can strictly be used only m cases where the shell is shallow. However, 
reasonably good accuracy is seen even when the elements are used to analyze 
shells that are only marginally shallow The user is cautioned, however, to be 
careful while interpreting results obtained when the shell analyzed is very 
deep Due to the excessive computation time associated with such calculations, 
the transverse shear flexibility is not taken into account m the element 
formulation. The element can be used m the statics, normal modes and 
differential stiffness rigid formats. 

Derivation of Element Properties 

Element geometry Rectangular Cartesian coordinates are used in the 
formulation An arbitrary triangular element is shown m figure A4 X , 

Y , and Z are the basic coordinates, x , y , and z are the local 
coordinates The grid points of the element are numbered m counterclockwise 
direction as shown m the figure. 

The lengths a , b , and c shown m figure A4 can be easily evaluated 
from the basic coordinates (Xi, Y x , Z x ) , (X 3 , Y 3 , Z 3 ) and (X 5, Y s , Z 5 ) of the 
vertices of the triangle 
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Displacement Field The u(x,y) and v(x,y) displacements are assumed to 
vary quadratically with position on the plane of the element, while displacement 
w(x,y) within the triangular element is assumed to vary as a qumtic polynomial 
m the local coordinates. 


u(x,y) = ai + a£X + a 3 y + a 4 x 2 + a 5 xy + a 6 y 2 


v(x,y) = a 7 + a 8 x + a 9 y + a 10 x 2 + a u xy + a 12 y 2 


w[x,y) = a 13 + a 14 x + a 15 y + a 16 x 2 + a 17 xy + a 18 y 2 

CD 

+ a ig x 3 + a 20 x 2 y + a 21 xy 2 + a 22 y 3 + a 24 x 3 y 
+ a 25 x 2 y 2 + a 26 xy 3 + a 27 y 4 + a 28 x 5 + a 29 xV 
+ a 30 x 3 y 2 + a 31 x 2 y 3 + a 32 xy 4 + a 33 y 5 


In concise form, u , v and w can be written as 


33 m n 

u = £ a i x 1 y 1 a i = m i = n i = Oi i = ? to 33 (2) 

i=l 


33 p q 

v = £ b x 1 y 1 b. “ P ! = q i = °, i " 1 t° 6 ( 3 ) 

i=l 

i = 13 to 33 


33 r s 

w = V cx 1 y 1 c = r = s = 0, i = 1 to 12 (4) 

*-i i 1 ill* 

1=1 


The detailed derivation of the stiffness matrix for the triangular shell 
element follows closely that for the TRIM6 and TRPLT1 elements. Hence, only 

the salient features of the derivation are given m this section 

✓ 

The geometry of the shell surface is approximated by a quadratic poly- 
nomial in base coordinates 
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z (x,y) = hi + h 2 x + h 3 y + h 4 x 2 + h 5 xy + h 6 y 2 


Hence the curvatures of the shell surface are 


z , = 2h 4 

9 xx H 


z > - hs 

xy ° 


J ’yy ■ 2h6 


The membrane thickness of the shell element is assumed to vary linearly 
over the surface of the element, i.e.. 


t. u 
1 .. 1 


- E 


The bending thickness of the shell element is also assumed a similar linear 


variation 


3 t* u! 

t K = £ d* x 1 y 1 


Following the shallow shell theory of Novozhilov (ref. 2) , the membrane 
strains m the shell are given by 


9u 

— z, w 

9x ’xx 


33 . 

= & K a * 


m -1 n. 
i l 
x y 


v s 


2h 4 c i x 1 y ^ 


3v 

- z , w 
9y yy 


P q.-l v s. . 

x 1 y 1 -2h6 c t x 1 y 1 ^ 
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3u 3v 

£ = "5 — + -r — - 2 z w 

xy 9y 9x > xy 

33 . m n -1 p -1 q 

“ E ( n x a ± x 1 y 1 + p x b ! x 1 y 1 Ci2) 

V s \ 

-2h 5 x 1 y 

In the absence of transverse shear effects, the bending strains are given by 


X x 


3 2 w 

3x 2 


33 v -2 s 

E v , (v - 1) c x 1 y 1 

i=l 1 1 


(13) 


X 


y 


9 2 w 

3x 2 


E s (s - l) c X 1 y 1 

i=l 1 1 1 


(14) 


2 d 2w 

X xy ** Z 3x9y 


33 

= E 

1=1 


2 v s. c 


v -1 s -1 

l l 

x y 


(15) 


Following the procedure outlined m sections 5.8.6 and 5.8.7, the jth 
column of the ith row of the generalized stiffness matrix is obtained as 

Ki 3 ° £ CGu (”i “j d k F(m i * m j * *k - 2 > n i * ", * “k J 

- htf. m. d k F(m x + + \ - 1, n ± + + u R ) 

- h 4 d k F (m^ + r x + t k - 1, n^ + s x + u R ) 

- hj ^ F(r 1 * Tj - t k , s x + Sj + up) 

( (16) 

+ d k F( P;l + + t k , q x + + U R - 2) 

- h 6 q x d k F( Pi + ^ + t k , q x + Sj + - 1) 

- h 6 q. d k FCr x + P 3 + V s x «■ qj + ^ - 1) 

+ ^6 d k F (r_^ + r^ + t k , S 1 + + u lP i ) (continued) 
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* n, d k F [ra^ + in + t k , * n + a, - 2) 

+ ”l Pj d k ^“i + Pj * \ - X > \ + ‘b + “k - 15 

' hs "i d k F ( m i + r j * l k’ n i * s j + u k ‘ W 

* P x n j \ FCPi ♦ ♦ t k - 1 , + n 3 + ^ - 1 ) 

+ P t Pj d k F(p. + p. - t k - 2, q L - qj + 

- h 5 Pi \ p (Pi + r 3 * \ - l " * s j * V 

- h 5 nj d k FC^ * m. + V »! + n, * ^ - U 

- hs P 3 d k F(r. + p. * \ - 1. s, ♦ Ij * u,J 

+ 4 \ FCr 1 * rj * t fc , Sj * Sj * u^) 

» On^ q d k F(m t ♦ P, * \ - 1. \ * ’ll * “]; - D 

- h 6 m i d k FCiik + r^. + t k 1, + 1^) 

" h b <k, d k PC^ + Pj + t k , s. + q. + u k - 1) 

+ 2h 4 h 6 FC^ + + t k , s ± + 5^ + v^) 

+ «i m 3 d k F( Pi + “j + \ - X > \ + n j + u k - 15 

- h 4 qi d k F( Pl + + t k , q x + s. * - 1 ) 

- h 6 m .j d k F(r. + + ^ - 1, + ru + 

* Gia^ d k F(m i + - 1, n. + - 1) 

+ m i p 3 d k FCm x + Pj + hi - 2 * \ + + V 

- hs m. d k F(m. + + t k - 1, m + s. + u R ) 

- h 4 d k FCr x + + t k , s_ L + + u k - 1) 

- h 4 Pj d k F(r. + Pj + t k - 1, s x + qj + u k ) 


C16) 


(continued) 
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+ 2hn h 5 d k F (r i + + t fc , -s i + 5^+1^) 

+ n x d k F + t R - 1, ^ - 1) 

- h„ n. d k F( Bl * r. * t k , n x * s, + \ - 1) 

* ?! m 3 d k Ff ?i * "j * h - 2. <1* + ", ♦ u k ) 

- h„ p x d k FCp 1 * ", + \ - 1. % * s , * <\) 

- h 5 in^ d k F[r. * mj ♦ ^ - 1, s t + nj + u*)) 

♦ G 2 3(q i n d k F( Pl + » t k , q x * n, * " k - 2 ) 

+ lx Pj d k F(p i + Pj * \ - !* % * <1, * “lc - D 

- h5 <>i d k FtF i + r , + V ’i * 5 j + "k - 2) 

- h 6 d k FCr x + + t k , s ± + + \ - 1) 

- h e P, \ FCr 1 * P, * \ - 1. s x * 1j + V 

+ 2h 5 h 6 d k F(r x + + t^, s i + + u fc ) 

+ n q d, F(m + p + t, , n + q + u, - 2) 

l n j k v l k’ l k 

- h 6 n x d k F(a 1 * r 3 * t fc , ^ * Sj * u k - 1) 

+ ?i 1, d k FC Pi + Pj * c k - l > % + + “k ' 1J 

- he Pl d k F( Pi tyVl, u k ) 

- h 5 d k F(r x * Pj ♦ t k , Si * q 3 t u k - 1))] 


* £ £ £ Cw d i, ,1 i, d v, < G » r x r , 

ki^i k^i k 3 =i 12 Kl * 3 1 J 


C" 


l)(r 3 - 1) 


F(r + r + t,' + t' + t' - 4, s + s 

i 3 k x k 2 k 3 1 J 


* "i, * u k 2 * V 


( 16 ) 


(continued) 
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A 


+ G 22 s s (s - l)(s. - 1) F(r + r + t.' + t,' + t, f , 

1 J 1 3 1 3 ki k 2 k 3 * 


+ s + s . + u,’ + u,' + u,* - 4) 

1 3 Tci TC2 K3 

+ ^4G 33 r i Tj s i + Gi 2 {r i (r x - 1) (s^ - 1) 

* r j s i Cr j - 1Ks i - 1)} j FCr i * r j * % * 'i 2 * % - 2 ’ 

" s i * s i + “ki + % * “k 3 ‘ 2) 

+ 2 Gi 3 {r i Sj (r. - 1 ) + ^ r ; s 4 - 1 )> FC^ . 

" Hu * % * % - 3 ’ s ! + s j * “k! + H z + u k 3 ' 15 

+ 2G 23 {r. s s (s - 1) + r. s- s. (s - 1)} F(r + r. 

j 1 j v 1 J iij'-j J < K x 3 

* C ki * % * % - 1 ' 3) 3 


( 16 ) 

(concluded) 


The generalized stiffness matrix can be transformed to the element and 
global coordinates by transformations similar to that for TRIM6 and TRPLT1 
elements . 


Equivalent Thermal Load Vector : 

The equivalent thermal load vector for the shallow shell triangular 
element consists of loads due to thermal expansion as well as due to thermal 
bending caused by variation of temperature with depth. The detailed derivation 
of the thermal load vector is similar to that used for TRIM6 and TRPLT1 

elements; hence, only the essential steps are given here. 

/ 

The vector of thermal strains is 



s xt 


a 

e l 


{e t } = ■ 

V 

> ss i 

a 

e 2 

( T - T ref) * <%> ^ " W 


‘ £ xyt ; 


a 

e 12 



(17) 


where {a g } = [u] -1 {a^} is a vector of thermal expansion coefficients, [u] 
is the strain transformation matrix given in equation (15) of page 5.8.4, {a^} 
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is the vector of thermal expansion coefficients m the material axis system, 

T - is the reference or stress-free temperature of the material and T is the 
temperature at any point (x,y) xn the element 

An applied stress vector which would produce the thermal strains is 


{£ t> - [ G J < E t> - C G el ( T * W 

The generalized equivalent thermal load vector 
{P gen } ■ TOT { U)t ^ 


(18) 



is obtained as 


(19) 


The strain vector {s} 


is given by 


(e) 


e 


3u 

-T z w - z 

x 

X 


3x xx 

A x 

e 



3v 

z w - z 

x 

y 


3y yy 

y 

. £ 


|H + |1. 2Z 

W - Z x 

xy 


3y 3x xy 

A xy 


( 20 ) 


where z , z and z are the curvatures of the shell surface and z is 
xx yy xy 

measured from the neutral surface of the plate. 

The temperature at any point (x,y,z), T , is given by 


T = T q + T' z (21) 

where T q is the mean temperature and T* is the thermal gradient. 

S 

'The following derivation to obtain the equivalent thermal load vector is 
given for the case of linear variation of thermal gradient over the planar 
coordinates, of the element; the values of the thermal gradient at the three 
vertices being defined as TJ , T^ and T£ . This capability is not 
operational m NASTRAN currently. The derivation, however, is valid for 
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cases with the same thermal gradient at the three vertices by setting 
and T£ equal to . Thus, T q and T’ o£ equation (21) vary over the 
element as follows: 


T Q = ei + e 2 x + e 3 y 


T* = e’ + e^x + e^y 


1 • C » y 


T > 
o 


3 


E 

1=1 



( 22 ) 

(23) 


(24) 


T* = 


3 

E 


1=1 



(25) 


The constants e^ , e 2 , e 3 and ej , e^ and e^ can be evaluated from 
the user supplied values of the mean temperature and temperature gradient at 
the vertices of the element; however, as stated earlier, only the capability 
of specifying a temperature gradient for the element is currently available 
and hence ej will be equal to the element temperature gradient and e£ 
and e| will be equal to zero. 


Substituting equations (10) through (15) into equation (20) and substituting 

for {e} and {cr, } in equation (19), the generalized equivalent thermal load 
£ ^ 

vector {P } is obtained as 
gen 
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33 m -1 n. v s 

l_, (m a x y - 2 h 4 c x y 
i=l 

v -2 s 

- z v (v - l) c x 1 y 

i i i 


/pt } _ _ j. 
iP gen^ " 9(a) 


33 Pi V 1 V i S i 

Z C q i b x x y - 2h6 c i x y 
i=l 

v s .-2 

- 2 Si (S 1 - 1) c i X 1 y 1 ) 

33 m n -1 P x -1 q 

V (n a x 1 y 1 +pbx y 1 
i=l 

v. s v -l s -1 

- 2h 5 c x X 1 y 1 - 2Z s x ^ X 1 y 1 )j 


[ 3 v. w v' w 1 "j 

^ C e x-*y^+e' x^y^z) dxdydz 

3=1 J ^ J 


Integrating over the thickness and noting that 


f f(Xjy) z dx dy dz = 0 
~t/2 


equation (26) reduces to 


{P t } 
gen 


9 {a} ff 


33 m -1 n vs 

E l 1 ou 11 

m i a i x y " 2h 4 c i x y 

i=l 


33 p q -1 Vs 

L % b , x 7 ~ 2h 6 c x y 

i=l 


33 » x n -1 Pi -1 q x 

Z-, n , % x y + P b X y 
i=l 


v s 

ii 

2h 5 x y 


V w 
x •* y ** 


**k \ 

d k x K y K ) dxdy 
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8{a} 


hff 


v . -2 s. i 

r 1 1 

V ! CV i - 15 C 1 * V 


v. s -2 

f T. 11 

s (s. - 1) c x y 
i v i i J 


v.-l s -1 

2v. s. c. x 1 y 1 
xii 


l%1 {« e } 


(28) 

(concluded) 


O o v! w'\ o o ^ r. +t. +t, 

C e: x 3 y j) v £ E ^AA, * 1 2 3 


3 3 3 


t, +t. +t. 


> 1 3 


k]=l k 2 =l k 3 =l 


ki“k 2 Tc 3 


• /VV^S dxdy 


The generalized equivalent thermal load vector will be obtained by performing 
the differentiation and integration operations of equation (28) and the final 
expression for {P^ en > will be similar to those obtained for the TRIM6 
and TRPLT1 elements, except that an additional expression involving the 
curvatures of the shell surface h^ , h 5 and hg will be added now. The 
generalized thermal load vector (Pg en ^ can then be transformed to the element 
and global coordinate system by the usual procedures. 


7.3.6. Differential Stiffness Matrix for Triangular Shell Element TRSHL 

The expression that is used for the energy of differential stiffness per 
unit area of the shell element consists of a part U£ due to out-of-plane 
motions and a part due to in-plane motions. The expressions for U£ 

and IP are the same as for plate elements and given in equations (18) and 
(19) of section 7.3.1.; the expressions for membrane strains will, however, 
involve the effects of coupling due to bending. Thus, 


U = IP + IP 
b m 


( 1 ) 


where 


U b 


t 

2 




— 9w 3w \ 
T xy 3x 3y / 


( 2 ) 
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and 


U' 

m 


-*{ 


0 (to 2 + 2to^ 


s ) 

xy' 


+ o 


K - 


2to 


£ ) 
xy 


+ 2t 


xy 


(e - 

y 


£ ) 
x 


(0 


(3) 


The stresses o , a and t at any point within the element is 
x y xy J r 

assumed to vary linearly, the values at the three corner grid points being 
used to evaluate the coefficients m the linear variation. 


<* x (x,y) = ei + e 2 x + e 3 y 

t 

° y (x,y) = fi + f 2 x + f 3 y 
^ (x,y) = gj + g 2 x + g 3 y 
In condensed form 


a 

x 


a 


y 







R 


x 


l 


S 


y 


i 


Also 


(0 

X 


3w 

3y 


to 


y 


3w 

3x 


(4) 

(5) 

( 6 ) 


(7) 


( 8 ) 


C9) 


CIO) 


( 11 ) 
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U) - 


1_ ( 3v 3w \ 
2 \3x ~ 3y/ 


( 12 ) 


e = 


3u 

Sx “ z ’ 


xx 


w 


(13) 


3v 

£ = -r z . w 

y ^y *yy 


(14) 


c = p. . !“ - 2z, w 
xy 3x 3y ’xy 


The thickness of the element t at any point 


(15) 


3 t v it 

t(x,y) = £ d, x K x 

i=l K 


(16) 


The jth column of the ith row of the generalized differential stiffness 
matrix is 


3 3 


K. = 


L y, fd, e„ r r Ffr + r + t, + R - 2, s + s. + u. + S„) 
k= 1 ^ k & i j ^ i j k H l J k. jr 

* d k f i s i s 3 F(r i W* s r 25 

* \ g i S 1 r j FCr i * r : * ^ + R >! - S 1 + s j * “k * s * - l) 

* d k % s 3 r i FCr i * * l k * R * - *• S 1 * s : - “k - s si - 15 

* 0.25 d k e s Pi Pj F( Pi ♦ P 3 * t k * \ - 2, q x + q., * ^ 

+ 0.25 d k e^ n i n^ + m^ + \ + \> \ * n j + \ + s jL ~ 2 ) 

- °- 25 d k e <t Pi n 3 F 'P! * “j * ‘k * R k - <k * n ] * “k * S « - 15 

- 0.25 d, e. n. p Ffm + p + t. + R. - 1, n + q. + u, + S. - 1) 

k£i r ]'-i^j k l ’ l n j k it 

+ ^k e A p i n j F ^ p i + ra j + t k + q i + n j + ^ + S £ ” (continued) 


(17) 
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+ d k e * ?! Pj F <Pi + p 3 + \ + h - 2 > \ + q 3 + u k + sj 
~ d k e £ Pi h 5 PCP X s £ ) 

" d k e * n x n j F(iu i + m 3 + \ + R v \ + n 3 + \ + s i ~ 2 ) 

- ^ n i P J F(ln i + Pj + \ + R £ - ^ n i + ^ + \ + S £ - 1) 

+ d k e Z \ h 5 F( \ + r j + \ + R £> n 2 + s 3 + \ + - 1) 

+ °‘ 25 d k f A Pi Pj F CPi + Pj + tfc + ** - 2, q x + q. * ^ + S £ ) 

" 0,25 d k £ * p i n 3 F( Pi + v t k + R r 1 ’ ( ii + y\ + y i) 

+ °* 25 d k f A n i n : p C ffi i + + \ * R £ , n x + ny u R + S £ - 2) 

" 0,25 d k f £ n ! Pj V( -\ + Pj + \ + R i - !> \ + q 3 - Ufc + S A - 1) 

- d k f * Pi n j F CP X + m y t k + \ - 1, q x + n 3 + u k + S A - 1) 

~ d k f £ p £ Pj F( Pi + Py t k + R £ - 2, q x + + u R + S £ ) U?) 

+ d k f * Pi h s + r y \ + R * - i. q t + ^ * u k + S £ ) 

+ d k f l n i n j F(m i + + h + V n ! + n j + \ + S £ - 2) 

+ d k f * n i P 3 p C\ + Pj + \ + R * - 1> \ ♦ q 3 + ^ - 1) 

' d k f A n i h 5 F (\ + T j + \ + R r n i + s j + \ + S £ - 1) 

+ 0,5 d k g £ \ Pj p CPi + py \ + \ - i, q x + q : + + S A - 1 ) 

+ ° 5 d k g £ % P i F ^i + Pj + \ + h ~ 1> \ + Rj + u k + S & ~ 1) 

" 0,5 d k g ^ \ n j F CPi + y y \ * n } + u k + s £ - 2 ) 

- 0-5 d k g £ n x F( P;j + a» 1 + t k + R £ , ^ + S £ - 2) 

°* 5 d k Pj hg F ^ V i + Pj + y + F jj, “ 1 j S i + + \ + (continued) 
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- 0.5 d k g t p. h 6 F(V. * P, - t k * Rjj - 1, 3j * q. . ^ ♦ S ( ) 

* °- 5 ^ g * "j h « F(I 1 U 

* °' 5 d k h \ h6 FCr 3 * "i + l k * R *’ s j * n i + \ * S jt - D 

- 0-5 d k Z t »i Pj FCm. ♦ Pj * \ ♦ Rj(. - 2, \ * q, * u*. ♦ sp 

- 0.5 d* g t TOj p t F(m. ♦ Pl ♦ t k + R t - 2, n + ^ ^ * S ) (17) 

(concluded) 

+ °* 5 d k % m i n 3 FCm i + m j + \ + \ ~ n x + n j + \ + S A - i) 

+ °* 5 ‘He g £ m 3 n i F(m i + *3 + \ + h ~ \ + n 3 + \ + S £ - U 

+ 0,5 ^ % P 3 hl * F(r i + Pj + \ + h " l > s i * <L, + \ + s £ ) 

+ °‘ 5 *k H Px ht > F < r j + Px + \ + R * ’ l > s j * <l x + \ + S t ) 

- 0.5 d k g z n. h^ F(r^ + m. + \ * \> \ * n. + + S £ - 1) 

- 0.5 d k g z n x h 4 FCr^ + m. + t k + R z , Sj + ^ ^ + S £ - 1)] 

The generalized differential stiffness matrix is transformed to the 
element coordinate system, basic coordinate system and the global coordinate 
system in the normal manner. 
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Section 15.2 


Modeling of Plate Structures Using TRPLT1 Elements 

The Figure 1, shown on page 15.2-3, is modeled using higher order 
triangular bending element, (Figure 2, page 15.2-3(a)), CTRPLT1 

Because of symmetry, the quarter section of the plate is discretized 
and detail of the discretization is given by the side of the modeled figure. 
Four different mesh sizes are used for each case. 

The central deflection is plotted in figures on pages 15.2-4 to 15.2-11 
and also given m Tables 1 and 2 on 15.2-3(b) and 15 2-3(c). 

Such high accuracy is obtainable for other plate structure problems 
using the TRPLT1 element m view of the use of the quintic displacement field 
for the displacement pattern in the element. 


15.2-2(a) (1/1/77) 
< 
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Table 1 Central deflection of simply supported rectangular plates 


b 

a 


= 2 . 


Number 

of 

Elements 

per 

Side 

N 

Concentrated 

Load at Center 

Uniformly Distributed Load 

Q-mesh 

P-mesh 

Q-mesh 

P-mesh 

2 

21.3344 

19.0681 

10.9158 

10.17408 

4 

17.7814 

17.1537 

10.1459 

10.0245 

8 

16.9230 

16.6984 

10.0776 

10.0548 

12 

16.7212 

16.6073 

10 1320 

10.1230 

Exact 

Solution 

16 

.5 

10. 

125 





























Number 

of 

Elements 

per 

Side 

N 

Concentrated 

Load at Center 

Uniformly Distributed Load 

Q-mesh 

P-mesh 

Q-mesh 

P-mesh 

2 

10.4294 

10.8878 

3.9168 

3.870672 

4 

8.4193 

8.0427 

2.7757 

2.7453 

8 

7.6242 

7.4392 

2.5791 

2.5738 

12 

7.4282 

7.3293 

2.5603 

2.5585 


i 



Exact 

Solution 


7.22 


2.54 

































Deflection Coefficient (a x io 3 ) 



Mesh Size (n) 


Central deflection of rectangular plate 
Case 2 (2-SS-U) 

15.2-5 (1/1/77) 



Deflection Coefficient (a x 10 3 ) 



Mesh Size (n) 

Central deflection of rectangular plate 
Case 3 (1-C-U) 

15 2-6 (1/1/77) 



Deflection Coefficient (a x 10 3 ) 



Mesh Size (n) 

Central deflection of rectangular plate 
Case 4 (2-C-U) 

15.2-7 (1/1/77} 




Mesh Size (n) 


Central deflection of rectangular plate 
Case 5 (1-SS-C) 

15 2-8 (1/1/77) 
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Mesh Size (n) 


Central deflection of rectangular plate. 
Case 6 (2-SS-C) 

15 2-9 (1/1/77} 
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Mesh Size (n) 


Central deflection of rectangular plate 
Case 7 (1-C-C) 

15 2-10 (1/1/77) 










Section 15.4 


Modeling Membrane Plate Using TRIM6 Element 

In the same figure 2 , the cantilever beam is discretized using the 
linear strain triangular membrane element TRIM6. Discretization and result 
of the corresponding displacement is shown on page 15.3-3. 

The cantilever beam shown on page 15.3-3 (figure 2) is divided into 
eight equal triangular (TRIM6) elements. The displacement pattern obtained 
using this mesh coincides with the exact one. Such high accuracy is 
obtainable for other membrane plate problems using the TRIM6 element m 
view of the quadratic displacement polynomial for the element. 


15.3-4 (7/1/76) 
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Nondimensional Deflection, v/c 


Modeling Errors in Membrane Plate Elements 
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©NASTRAN TRIM6 
QDMEM1 
QDMEM 



Figure 2. Deflection of cantilever beam idealized by QDMEM1 and TRIM6 elements. 

15 3-3 (7/1/76) 


84 



APPENDIX B 

Updates to the NASTRAN Users' Manual 
for the addition of TRIM6, TRPLT1 and TRSHL elements 
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STRUCTURAL ELEMENTS 


1.3.5. Plate Elements 

NASTRAN includes two different shapes of plate elements (triangular 
and quadrilateral) and two different stress systems (membrane and bending) 
which are uncoupled. There are in all a total of thirteen different forms of 
plate elements that are defined by connection cards as follows 

1. CTRMEM - triangular element with finite m-plane stiffness and zero 
bending stiffness. 

2. CTRIM6 - a triangular element with finite m-plane stiffness and zero 
bending stiffness. 

3. CTRBSC - basic unit from which the bending properties of the other 
plate elements are formed. 

4 CTRPLT - triangular element with zero m-plane stiffness and finite 
bending stiffness. 

5. CTRPLT1 - higher order bending element — a triangular element with 
zero m-plane stiffness and finite bending stiffness. 

6. CTRIA1 - triangular element with both m-plane and bending stiffness. 

It is designed for sandwich plates which can have different materials referenced 
for membrane, bending and transverse shear properties 

7. CTRIA2 - triangular element with both m-plane and bending stiffness 
that assumes a solid homogeneous cross section. 

8. CQDMEM - quadrilateral element consisting of four overlapping CTRMEM 
elements. 

9. CQDMEM1 - an isoparametric quadrilateral membrane element. 

10 CQDMEM2 - a quadrilateral membrane element consisting of four 
nonoverlappmg CTRMEM elements. 

11. CQDPLT - quadrilateral element with zero m-plane stiffness and finite 
bending stiffness 


1.3-5 (4/1/73) 


S, 


\ 

PA®® 
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12. CQUAD1 - quadrilateral element with both m-plane and bending stiff- 
ness. It is designed for sandwich plates which can have different materials 
referenced for membrane, bending and transverse shear properties. 

13. CQUAD2 - quadrilateral element with both in-plane and bending stiff- 
ness that assumes a solid homogeneous cross section. 

Theoretical aspects of the plate elements are treated in Section 5.8 of the 
Theoretical Manual. 

In addition, a shallow shell element, CTRSHL is also available. The 
elements and the coordinate systems are shown m figures 14(a), (b) and (c) 


1.3-6 (3/1/76) 
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BULK DATA DECK 


Input Data Card CTRIM6 Triangular Element Connection 

Description : Defines a linear strain triangular membrane element (TRIM6) 

of the structural model. 

Format and Example : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

CTRIM6 

EID 

PID 

Gl 

G2 

G3 

G4 

GS 

G6 

A1 

CTRIM6 

220 

666 

100 

110 

120 

210 

220 

320 

+C2 


+BC 

TH 









+C22 

9.0 










Field 


Contents 


EID 

PID 

GI thru G6 


TH 


Element identification number (integer > 0) . 

Property identification number (integer > 0) . 

Grid point identification numbers of connected 
points (integers > 0; Gl f G2 f G3 f G4 f G5 i G6) . 

Material property orientation angle in degrees 
(Real) . The sketch below gives sign convention 
for TH. 
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Remarks : 


The grid points must be listed consecutively going around the 
perimeter in an anticlockwise direction and starting at a vertex. 


2. Material properties (if MA.T2) and stresses are given m the 

(x . y ) coordinate system shown m the sketch. 
hi 'm 

3. G2, G4, and G6 are assumed to lie at the midpoints of the sides 
The locations of these grid points (on GRID Bulk Data cards) are 
used only for global coordinate system definition, GPWG (weight 
generator module), centrifugal forces, and deformed structure 
plotting . 

4. Continuation card must be present. 

5. Element identification numbers must be unique with respect to all 
other element identification numbers 
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BULK DATA DECK 


Input Data Card PTRIM6 Linear Strain Triangular Element Property 

Description : Defines the properties of a linear strain triangular membrane 

element TRIM6 . 

Format and Example : 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

PTRIM6 

PID 

MID 

Tl 

T3 

T5 

NSM 




PTRIM6 

666 

999 

1.17 

2.52 

3.84 

8.3 





Field 

PID 

MID 

Tl, T3, T5 
NSM 


Contents 

Property identification number (integer > 0} . 
Material identification number (integer > 0) . 
Thickness at the vertices of the element (Real) . 
Nonstructural mass per unit area (Real) . 


Remarks : 

1. For structural problems, the material may be MAT1 or MAT2. 

2. The thickness varies linearly over the triangle. If T3 or T5 
is specified 0.0 or blank, it will be set equal to Tl. 

3. All PTRIM6 cards must have unique property identification numbers. 
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BULK DATA DECK 


Input Data Card CTRPLT1 Triangular Element Connection 

Description Defines a triangular bending element (TRPLT1) of the structural 
model. 

Format and Example 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

CTRPLT 

EID 

PID 

G1 

G2 

G3 

G4 

G5 

G6 

+abc 

CTRPLT 

160 

20 

120 



40 

70 

110 

+ABC 



TH 










16 2 










Field Contents 

EID Element identification number (integer > 0) 

PID Identification number of a PTRPLT property card 

(Default is EID) (integer > 0) . 

Gl, G2, G3, G4, G5, G6 Grid point identification numbers of connection 

points (integer >0* G1 / G2 1 G3 f G4 i G5 i G6) . 

TH Material property orientation angle m degrees (Real) 

The sketch below gives the sign convention for TH. 


Y 
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Remarks * 


1. Element identification numbers must be unique with respect to all 
other element identification numbers. 

2. Interior angles must be less than 180°. 

3. The grid points must be listed consecutively going around the 
perimeter m an anticlockwise direction and starting at a vertex. 

4. Continuation card must be present. 
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BULK DATA DECK 


Input Data Card PTRPLT1 Triangular Plate Property- 

Description* Used to define the bending properties of a triangular plate 
element. Referenced by CTRPLT1 card. No membrane properties are included. 

Format and Example 


1 

2 

3 

4 

5 

6 

7 

3 

9 

10 

PTRPLT1 

PID 

MIDI 

11 

13 

15 

MID 2 

TS1 

TS3 

+abc 

PTRPLT1 

15 

25 

20 

30 

40 

35 

3.0 

1.15 

+PQR 


+abc 



Zll 



Z23 

Z15 

Z25 


+PQR 

1.0 

-1 0 

1 5 

-1.5 

2.0 

-2.0 

2.5 

-2.5 



Field 


Contents 


PID 


Property identification number (integer > 0). 


MIDI 


Material identification number for bending (integer > 0) 


II, 13, 15 


MID2 


Area moment of inertia of the element per unit width at 
the vertices 1, 3, 5 of the element (Real > 0.0) 



where Tj , T 3 , T 5 are the thickness of the element 
at the vertices 1, 3, 5. 

Material identification number for transverse shear 
(integer > 0) . 


TS1 , TS3, TS5 


Transverse shear thickness (Real > 0 0) at the 
vertices 1, 3, 5 of the element. 


NSM 


Nonstructural mass per unit area (Real) 


Zll, Z21, Z13, Z23, Fiber distances for stress computation at grid points 

Z15, Z25 Gl, G3, G5, respectively, positive according to the 

right-hand sequence defined on the CTRPLT1 card (Real) 
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Remarks : 


1. All PTRPLT1 cards must have unique property identification numbers. 

2. If TS1 is zero, the element is assumed to be rigid in transverse 
shear. 

3. If TS3 or TS5 is 0.0 or blank, it will be set equal to TS1. 

4. If 13 or IS is 0.0 or blank, it will be set equal to II. 

5. The stresses at the centroid will be computed at the top and bottom 
fibers . 




A 



BULK DATA DECK 


Input Data Card CTRSHL Triangular Shell Element Connection 

Description : Defines a triangular thin shallow shell element CTRSHL) of 

the structural model. 

Format and Example 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

CTRSHL 

EID 

PID 

Gl 

G2 

G3 

G4 

G5 

G6 

+abc 

CTRPLT 

160 

20 

120 

10 

30 

40 

70 

110 

+ABC 


+abc 

TH 









+ABC 

16.2 










Field 


Contents 


EID 


Element identification number (Integer > 0) 


PID 

Gl, G2, G3 
G4, G5, G6 


Identification number of PTRSHL property card 
(Default is EID) (Integer > 0) 

Grid point identification numbers of connection 
points (Integer >0: Gl $ G2 f G3 £ G4 i G5 £ G6) 


Material property orientation angle in degrees (Real) - 
The sketch below gives the sign convention for TH. 



Remarks 

1. Element identification numbers must be unique with respect to all 
other element identification numbers 

2. Interior angles must be less than 180°. 
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3. The grid points must be listed consecutively going around the 
perimeter in an anticlockwise direction and starting at a vertex. 

4. Continuation card must be present. 
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BULK DATA DECK 


Input Data Card PTRSHL Triangular Shell Property 

Description Used to define the bending properties of a triangular shell 
element. Referenced by the CTRSHL card. 

Format and Example 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

PTRSHL 

PID 

MIDI 

Tl 

T3 

T5 

MID2 

11 

13 

1 

PTRSHL 

10 

20 

3.0 

6.0 

4.0 

30 

2.25 

18.0 

+PQR 


+abc 

15 

■sa 

TS1 

TS3 

TS5 

NSM 

Zll 

Z21 

+def 

+PQR 

5.33 

40 

2.5 

5.0 

3.5 

50 


-1.5 

+STU 


+def 

Z13 

Z23 

Z15 

Z25 





i | 

+STU 

3.0 

-3.0 | 

2.0 

-2.0 







Field 


Content 


PID 

MID 


Tl, T3, TS 
MID2 


II, 13, 15 


Property Identification number (Integer > 0) . 

Material identification number for membrane effect 
(Integer > 0) . 

Thickness for membrane action at vertices 1, 3, 5 of 
the elements (Real > 0.0). 

Material identification number for bending effects 
(Integer > 0) . 

Area moments of inertia of the element at the vertices 
1, 3, 5 of the element (Real > 0.0) 


MID3 

TS1, TS3, TS5 


Material identification number for transverse shear 
(Integer > 0) . 

Transverse shear thickness (Real > 0.0) at the vertices 
1, 3, 5 of the element. 


NSM 


Non- structural mass per unit area (Real) . 


Zll, Z12, Z13, Fiber distances for stress computation at grid points 

Z23, Z15, Z25 Gl, G3, G5, respectively, positive according to the right- 

hand sequence defined on the CTRSHL card (Real > 0.0). 
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Remarks : 


1. All PTRSHL cards must have unique property identification numbers. 

2. If T3 or T5 equal to 0.0, or blank, they will be set equal to Tl. 

3. If 13 or 15 equal to 0.0, or blank, they will be set equal to II. 

4. If TS3 or TS5 equal to 0.0, or blank, they will be set equal to TS1. 

5. If TS1 is 0*0, or blank, the element is assumed to be rigid in 
transverse shear. 

6. The stresses at the centroid will be computed at the top and bottom 
fibers . 
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APPENDIX C 


Updates to the NASTRAN Programmer 1 s Manual 
for the addition of TRIM6, TRPLT1 and TRSHL elements 



4.87.21. TRIM6: Linear Strain Triangular Element 


4.87.21.1 Input Data for TRIM6 Element 


EST entries for TRIM6 are 


Symbol 


Description 


Element Identification Number 


SIL 1} SIL 2j 


Scalar indices of connected grid points 
Anisotropic material orientation angle 


Mat ID 


Material Identification Number 


Tl, T3, TS 


Thickness of comer grid points 


Nonstructural mass per unit area 
Local coordinate system numbers and 


1 = 1, 6 


location coordinates m the basic 


system for the connected grid 


Z^ j points 

T01, T02, T03, T04, TOS, T06 Temperatures at the grid points 

2. Coordinate system data 

The numbers N , X. , Y and Z are used to calculate 3 by 3 

ill i 

basic-to-global coordinate transformation matrices [T J for points 
i = 1, 2, 3, 4, 5 and 6 . 

3 . Material data 


Symbol 


Description 


3x3 stress -strain matrix 


Mass density 


a , a , a 
x y xy 








Thermal expansion coefficients 


erence temperature 
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g e Structural damping coefficient 

cr^, 0 c> a s Stress limits for tension, compression 

and shear 

4.87.21.2 Basic Equations for TR1M6 

1. The element coordinate system is defined by the following equations 


( x 3 - X X \ 

Y 3 - Yi 1 
z 3 “ Z 1 / 

( X 5 - Xj \ 
{Vi 5 > = < Y 5 - Y X i 

l z 5 “ Z 1 J 


CD 


(2) 


{Vi 3 } 

{i} = iwnn 


(3) 


{i} x {V 13 } 
{k} = |{i} x {V 13 }| 


(4) 


ij) = {k} x U} 


( 5 ) 


2. The displacement transformation matrix from basic coordinates to 
m-plane coordinates is : 



( 6 ) 


3. The local (element) coordinate system of the element is as follows. 
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The x-axis is obtained by joining grid points 1 and 3 of the element. 

The y-axis is the perpendicular from grid point 5 to the x-axis (line 
j ommg grid points 1 and 3) . 

Depending upon the location of grid point 5 relative to grid points 1 
and 3, 3 cases of triangle orientation are possible: (refer to fig. 4.87.21.1) 

Case I Acute angles at grid points 1 and 3 


c = |{ 1 } X {V 15 } | 


(7) 

b = {i} • {V 15 } 


(8) 

a = j{V 13 >| - b 


C9) 

Coordinates of points are 

u a - b 

a 


X! = -b ; x 2 = — 2 — * x 3 = a , x 4 = 

2 > *5 “ 0 , 

(10) 

b 

X 6 = - j 

n = o , yz = o ; y 3 = o , n = f , 

75 = C , y 6 = j 

(ID 

Case II. Obtuse angle at grid point 3 

c = |{i} x {Vi 5 }| 


(12) 

b = {i} • {V l5 } 


(13) 

a = b - | (V 13 }| 


(14) 


Coordinates of points are 

, -a - b 

X 1 - -b , x 2 = — 2 » x 3 = ~ a > x 4- 

b 

x 6 = - 2 
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12 3 12 3 

C-b,0) C-a,0) (b,0) (a,0) 

Case II. Obtuse angle at grid point 3. Case III. Obtuse angle at grid point 3 


Figure 1 Triangular element shapes. 
4.87-21.1 


1Q8 


( 16 ) 


yi * 0 > /2 = o ; y3 = ° ; y>* = §■ > ys = c , Y6 = j 

Case III Obtuse angle at grid point 1 


= [{l} X {Vl5>| 

( 17 ) 

= {i> • (Vis) 

( 18 ) 

= j ‘CV 1 3> | + b 

( 19 ) 


Coordinates of points are 




a + 

b 



a 


b 


xi = b 3 

■ x 2 

2 

9 

x 3 = 

a , x 4 

" 2 3 

x 5 = 0 ; 

x 6 = 2 


yi = o j 

: 72 

= 0 ; 

ys 

= 0 ; 

CM 

O 

II 

X 

> Ys 

= c ; ye 

c 

" 2 


4. The matrix [_Hj relating grid point displacements and the generalized 

coordinates (in the equation (u) 

= H 

{a}) is 

given by 




"Hi 

| 0 








Lh] = 

— 

J 







(20) 


0 

\ Hi. 








where 











1 

*1 

yi 

x l 

x iyi 

yf" 





1 

x 2 

y 2 

x 2 

x 2 y 2 

yi 




[Hi] = 

1 

x 3 

y3 

2 

x 3 

x 3y3 

2 

y3 







2 


2 



(21) 


1 

x 4 

y* 


x 4 y 4 

Yh 





1 

x 5 

ys 

2 

x 5 

x 5y5 

2 , 
ys 








2 


2 





1 

x 6 

ye 

x 6 

x eye 

ye 
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5. The matrix [b] relating strain vector to the generalized coordinates 
(m the equation {e} = '[b] {a}) is given by 


B = 


010 2xy0000000 

000*0000010x 2y 
00 1 0x2y0 1 0 2xy0 


( 22 ) 


4.87.21.3 Stiffness Matrix Calculation for TRIM6 (Subroutine KTRM6S and KTRM6D) 

The polynomial expressions for variation of u , v and t within the 
element are 


12 


m. n 
i i 


u = Z a , x 1 y 

i=l 


(23) 


12 p q 

V 1 u 11 

= Z b i x y 

i=l 1 


(24) 


r, s. 
i i 


t = Z c 4 x 1 y 
i=l 


(25) 


The values of m , n , p , q. ,r,s are 

i ’ i ' r i * n i 5 i l 


mi = 0 ; m 2 = 1 , m 3 = 0 ; m^ = 2 ; m 5 = 1 , m 6 = 0 , 


my to mi 2 = 0 


(26) 


n i “ 0 ; n 2 = 0 ; n 3 = 1 ; n 4 = 0 ; n 5 = 1 , n 6 = 2 , 
n 7 to n i2 = 0 


(27) 


Pi to P6 = 0 , p 7 = 0 , p 8 s 1 , P9 = 0 , p 10 = 2 , 

Pll = 1 * Pl2 " 0 


(28) 
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qi to q 6 = 0 ; q 7 = 0 , q 8 = 0 ; q 9 = 1 ; qio = 0 


( 29 ) 


qu 

= i > qi2 - 2 


ri = 0 ; 

r 2 = 1 , r 3 = 0 

(30) 

si = 0 ; 

s 2 = 0 , s 3 = 1 

(31) 

P 

rt 

O 

a 12 = 0 ; bi to b 6 = 0 

(32) 


The coefficients a^ to a 6 and b 7 to b 12 are generalized 

coordinates of the element and can be evaluated once the displacement vector 

is known. 

The coefficients ^ , c 2 and c 3 can be evaluated from the specified 
thicknesses tj , t 3 and t5 of the 3 comer grid points and the geometric 

dimensions a , b and c of the element 


t*a + t 3 b 
Cl (a + b) 


(33) 


t3 ~ *1 
Cz “ (a + b) 


(34) 



(35) 


The elements of the symmetric portion of the stress -strain matrix 
[G e ] are denoted by G n , G 12 , G 13 ^ G 22 , G 23 , G 33 . 

A formula for the integral of the type x ra y n taken over the area of 
the element is 


// 


m n , , .. n+1 , m+1 , .m+1. 

x y dx dy = F(m,n) = c (a - (-b) } 


min? 

(m + n + 2) * 


(36) 
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The equation used in the stiffness matrix generation in generalized 
coordinates is 


(k 


) = T c, [Giim m F(m + m + r, - 2, n +n + s, ) 

lj-’gen k L 11 1 ] v i j k 5 l j k' 

+ G 22 q i q :J F (p i + Pj + r k , q x + + s k - 2) 

+ G 33 {n i n F (m^ + + r k , + n^ + s k - 2) 

* Vj F + r k - 2 > - 1 - ^ * s k» 

+ (6330^ + ^2^) F Cm- + P 3 + r k - i, n i + q^ + s k - 1) 

+ (233^.^ + G l2 m.q.) F (nu + Pj. + r k " 1* n -, + q i + s k “ 13 

+ G ! 3 { Cm^ n x + nun^) F (n^ + m. + r k - 1, n i + n^ + s^ - 1) 

+ m j P i F (m 3 + p i + r k ~ 2 ’ n 3 + q i + V 

+ F (m x + Pj + x k - 2, n i + q^ + s^)} 

+ G 23 {(P i q j + P j q i ) F (p L + + r k - 1, q x + q^ + s k - 1) 

- Vj F ( ”i * Pj + r k’ n i - 1, * s k - 2J 

* n j q i F (m. * Pl ♦ r k , r ; * * s* - 2)}] 

The stiffness matrix in global coordinates is 


( 37 ) 


M - [B] W T [h _ 1 ] T [kl [h' 1 ] [T] [Ef 


gen 


( 33 ) 


For use m the overall structural matrix, the 3x3 fc partition 
of the stiffness matrix [k] corresponding to grid point i and connection 
point j is expanded to 6 x 6 to form 




(39) 
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4.87.21.4 Mass Matrix Calculation for the TRIM6 Element (calculated m the 
stiffness subroutine KTRM6S and KTRM6D) 

The mass is generated by the fol lcwing algorithm 


( X 3 - Xj ' 
Y 3 - Y X > 
Z 3 - Zj , 

1 X 5 - Xi \ 
Ys - \ 

Z 5 - Zi J 


(40) 


(41) 


The area is 


A = j | {V 13 > x {V 15 }[ 


(42) 


Volume 


V = ci F (0,0) + c 2 F(1,0) + c 3 F (0,1) (43) 

where Ci , c 2 , c 3 [see eq. (33), (34), (35)] are the constants in the 
thickness equation of the element [eq. (25)J and zero factorial has a value 
of 1. The mass at each point is 

m = -i- (pV + Au) (44) 


which is -g- the total mass. 

For each point the diagonal mass matrix m element coordinate system 
at all the grid points is 
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so that [u ] the element mass matrxx has [jn. 3 matrices arranged 

0© X 

diagonally. 

The mass matrix in global coordinate system is obtained as 

["„] = fe] [t] T [M e „] [t] [b] T (46) 

4.87.21.5 Element Load Calculations for TR1M6 (Subroutine TLi5DM6) 

The temperature within the element is assumed to vary bilmearly 

3 r. s 

T = £ d * V 1 (47) 

i=l 

with 

ri = 0 ; r 2 = 1 , r 3 = 0 (48) 

and 

Si = 0 , s 2 = 0 and s 3 = 1 (49) 

The coefficients dj , d 2 and d 3 are evaluated from the specified 
temperatures T 01 , T 03 , T os at the three corner grid points (obtained 
from the GPTT data block) and the reference temperature Tq of the element 
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Toia + To3b 


dl (a + b) 

(50) 

T 03 “ T01 


d 2 (a + b) 

(51) 

d 3 = “■ [t<) 5 - d l] 

(52) 


The constant dji is modified by the reference temperature, T 0 , 

di = di - Tq . The ith element of the generalized load vector {P } is 

gen 

3 3 

{ Vgen = gj S °k d 2 t G u m l F K + r k + \ - l * \ + s k * V 

+ G 22 F (p i + r k + t^, q x + + u £ - 1) 

+ {n i F (m i + r k + V n i + s k + u £ - 


(53) 


+ p i F (p i + r k + h ~ l > *1 + s k + 


where 


G 11 = G 11 cti + G 12 «2 + G 1 3 «12 

G 22 = G 12 + G 2 2 <*2 + G 2 3 a 12. 

G 33 s G 13 a l + g 23 <*2 + g 33 «12 


The generalized equivalent load vector {P } is transformed to 

gen 

load vector {P } m local element coordinates and to load vector {P } 
e g 

m global grid-point coordinates by the following transformations 


{P } = [h->] t {P } 


(54) 
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{P g } - [E] [t] t {P e } 


(55) 


{P } is a 18 x 1 vector, 
g 

The forces axe placed in the PG load vector data block. 

4.87.21.6 Element Stress Calculations for TRIM6 Element (Subroutine STRM61 
and STRM62 of module SDR2) 

1. The relationship between strain and generalized coefficients is 

{e} = [b] {a} (56) 

where 

010 2x^0000000 

[b]= 000000001 Ox 2y 

.00l0x2y0102xy0 

The transformation from displacements to stress is* 

iV] = [g 6 ] [b] [h- 1 ] [e] t [t] 

The temperature to stress relation is 

{s t > = “[ G e ] {ct} ( S9 ) 

where 
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for isotropic materials, {a} is input by the user for anisotropic materials 
and corrected for material angle by 


a = 



(61) 


2. Calculations performed by STRM62 (Phase 2 calculations) 
The equation for stress is 


a 

x 


o 

y 

a 

xy 


1 



* <s t > ct 3 


T ) 
o' 


(62) 


where is the loading temperature for the point where stress is evaluated 

(3 comer grid points and centroid) and is obtained from the GPTT data block. 
The temperature of the centroid is taken as the average of the grid point 
temperatures . 

The principal stresses are 


Oj = 


a 2 = 




0 = 


2 arctan 



m degrees 


where 0 is limited to: -90° < 0 < 90° 

The maximum shear is 


(63) 


(64) 


(65) 


T 



( 66 ) 


117 



The stresses are output for 4 points for every element. 3 comer grid 
points and the centroid. 
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4.87.22. TRPLT1 Higher Order Plate-bending Element 


4.87.22.1 Input Data for TRPLT1 Element 


1 


EST entries for TRIBI are. 
Symbol 


E1D 

SILi, SIL2, . . • , SILg 
0 

Mat ID^ 

Mat ID5 
II, 13, IS 


TS1, TS3, TS5 


Zll, Z21, Z13, Z23, 
Z1S, Z25 
N 


X 


1 - 1 , 6 


1 / 

TEMP 


Description 

Element Identification Number 

Scalar indices of connected grid points 

Anisotropic material orientation angle 

Material Identification Number for bending 

Material Identification Number for shear 

Area moment of inertia per unit width at comer 

*3 *3 t 3' 

t- 1 c 3 15 

grid points II = , 13 = , 15 = jj 

Effective thickness for transverse shear 

at corner grid points 

Nonstructural mass per unit area 

Distances Z1 and Z2 for stress 

calculation at 3 comer points 

Local coordinate system numbers 

and location coordinates in the 

basic system for the connected 

grid points 

Element temperature 


2. Coordinate system data 

The numbers N. , X , and Z are used to calculate the 3 by 3 
1 3 1 1 J 

basic-to-global coordinate transformation matrices []T ] for points 
i=l, 2, 3, 4, 5, 6, (via subroutines TRANSD or TRANSS) . 
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3. Material data 


Symbol 


&] 


For mat. 


ID, 


a , a , a 
x y xy 

TO 




i 0 * a , a 
\ t’ c J s 


For mat. ( G 


Description 
3x3 stress -strain matrix 
Mass density 

Thermal expansion coefficients 

Reference temperature 

Structural damping coefficient 

Stress limits for tension, compression 

and shear 

Shear coefficient 


ID. 


4.87.22.2 Basic Equation for TRPLT1 


1. The element coordinate system is defined by the following equations- 


{V 13 > = < 



( X 5 - Xi \ 


{Vi 5 } = 




Y 5 - *1 

v z 5 - z a ; 




A = i | {v 13 } x {v 15 }| 


CD 


C2) 


C3) 


^13* 

{l} = Uvi 3 >| 


( 4 ) 
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{l} X {V13} 
|{l} x {V^} 


( 5 ) 


{j} = {k} X {i> 


( 6 ) 


2 . The displacement transformation matrix from basic coordinates to 
in-plane coordinates is: 




k l 

k 2 

k 3 

0 

0 

0 


0 

0 

0 

M 

*2 

13 


.0 

0 

0 

3 l 

3 2 

33 . 


( 7 ) 


3 . The local (element) coordinate system of the element is as follows 

The x-axis is obtained by joining grid points 1 and 3 of the element. 

The y-axis is the perpendicular from grid point 5 to the x-axis (line 
joining grid points 1 and 3 ) . 

Depending upon the location of grid point 5 relative to grid points 1 
and 3 , three cases of triangle orientation are possible: (refer to fig. 4 . 87 . 21 . 1 ) 


Case I Acute angles at grid points 1 and 3 . 


c = |{i> x {V 15 >| 

(8) 

b = {1} • {V25} 

( 9 ) 

a = [ {V 13 >I - b 

(10) 

Coordinates of points are 



xj = -b 



*3 = a , 



x 5 = 0 , 


x 6 “ - 


b 

2 


( 11 ) 
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yi = o ; yz ~ o ; y$ = o ; y 4 = f- ; ys ■ c ; y 6 = f (12) 

Case II : Obtuse angle at grid point 3: 

c — |{i} x {V 15 }| (13) 

b = {1} • {V 15 } (14) 

a = b - | {V 13 }| (15) 

Coordinates of points are 

1 a ” b a a a 

xi = -b ; x 2 = § — ’ x 3 = 0 , x 4 = - j , x 5 = 0 , 

(16) 

b 

x 6 = “ 2 

yi 9 0 i yz =' 0 ; y3 s 0 ; y^ = § > ys = c ; y 6 = § (17) 

Case III * Obtuse angle at grid point 1. 

c = j{i) x { V 15)! (18) 

b = {i> • (V 15 > (19) 

a * |(V 13 }| + b (20) 

Coordinates of points are: 

, a + b a 

xi = b , x 2 = — ^ — ; x 3 = a , - x 4 = j ; x 5 = 0 , 

C21) 

b 

x 6 - 2 

yi “ 0 ; y 2 - 0 , y 3 = 0 ; y 4 = §- ; y 5 = c , y 6 = f (22) 
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The matrix [H] (for plates infinitely rigid in transverse shear) 
relating grid point displacements and the generalized coordinates (in the 
equation {u} - [h] {a}) is given by the matrix M . on the following 
page 

5 . The matrix M relating curvatures (for plates infinitely 
rigid m transverse shear) to the generalized coordinates in the equation 
fxi> = [b 2 ] {a}) is given by 

’o 0 0 2 0 0 6x 2y 0 0 12x 2 bxy 2y 2 0 0 20x 3 6xy 2 2y 3 0 0 

[b 2 ] = 000002 0 0 2xby 0 0 2x 2 bxy 12y 2 0 2x 3 6x 2 y 12xy 2 20y 3 

_0 0 0 0 2 0 0 4x 4y 0 0 6x 2 8xy 6 y 2 0 0 12x 2 y 12xy 2 8y 3 0 

(24) 

6. The matrix [bJ relating transverse shear strains {y} to the 
generalized coordinates (m the equation {y} = [Bj 3 {a» is a 2 x 20 matrix 
whose nonzero elements are as follows: 


Bid, 7) = 6A n (2Sa) 

Bi(l, 8) = 2 A 3 i (25b) 

Bi(l,9) = 2A 32 (25c) 

Bi( 1, 10) = 6A 1S (2Sd) 

Bj (1,11) = 24A u x (25e) 

Bi(l,12) = 6 (A 31 x + Any) (25f) 

Bi(l,13) = 4 (A 32 x + A 31 y) (25g) 

Bj (1,14) = 6(A 15 x + A 32 y) (25h) 

Bi( 1,15) = 24A 15 y (2Si) 
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0 

0 

i 

0 
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0 
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„ 2 
3y2 

0 
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. 3 
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0 
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2 

-y2 

0 

. 3 
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-3xfy 2 

-2x 2 y| 

3 

-y 2 

0 

c 4 

-5X2 

2 2 
-3x 2 y 2 

-2 x 2 y| 

4 

-y 2 

0 

K p 

L_> tarj 





2 


2 

3 

2 

2 

3 

4 

3 

2 2 

2 

4 

5 

3 2 

2 3 

4 

5 

a* 


1 

x 3 

y 3 

x 3 

x 3y3 

ys 

x 3 

x 3y3 

x 3y3 

y3 

x 3 

x 3y3 

x 3y3 

x 3y3 

y3 

x 3 

x 3y3 

x 3y3 

x 3y3 

y3 


0 

0 

i 

0 

x 3 

2y3 

0 

2 

x 3 

2 x 3 y3 

2 

syt 

0 

3 

x 3 

2x iy 3 

3x 3 y 3 

, 3 
4y3 

>0 

2x 3 y 3 

7 2 2 

3x 3 y 3 

4x 3 y 3 5 y 3 


[H] = 

0 

-1 

0 

- 2 x 3 

-y3 

0 

7 2 

-3x 3 

- 2 x 3 y 3 

2 

-y3 

0 

4 3 

-4x 3 

„ 2 

-3x 3 y 3 

- 2 x 3 y 3 

3 

-y3 

0 

-5x 3 

2 2 
-3x 3 y 3 

- 2 x 3 y| 

-yi 

0 




Y4 

2 


2 

3 

2 

2 

3 

4 

3 

2 2 

3 

4 

5 

3 2 

2 3 

4 

5 



1 

X 4 

x 4 

X4y4 

y4 

x 4 

x 4 y 4 

x 4 y 4 

y4 

x 4 

x 4y4 

x 4 y 4 

x 4 y 4 

y4 

x 4 

x 4 y 4 

x 4y 4 

x 4 y 4 

y4 



0 

0 

1 

0 

X4 

2 y 4 

0 

2 

x 4 

2 x 4 y 4 

syi 

0 

3 

x 4 

2 x 4 y 4 

3x 4 y 4 

. 3 
4x 4 

0 

2 x 4 y 4 

7 2 2 

3x 4 y 4 

4x 4 y 4 5y 4 



0 

1 

0 

- 2 x 4 

-y 4 

0 

-3xJ 

- 2 x 4 y 4 

-y 4 

0 

-4x 4 

-3x 4 y 4 

- 2 x 4 yJ 

-yj 

0 

-5x 4 

„ 2 2 
-3x 4 y 4 

-2x 4 y4 

-y 4 

0 




x 5 

ys 

2 


2 

3 

2 

2 

3 

4 

3 

2 2 

3 

4 

5 

3 2 

2 3 

4 

5 



1 

x 5 

x sys 

ys 

x 5 

xsys 

xsys 

ys 

x 5 

xsys 

xsys 

xsys 

ys 

x 5 

xsys 

xsys 

xsys 

ys 


i 

0 

0 

i 

0 

x 5 

2y5 

0 

x 5 

2X5y5 

3y| 

0 

x 5 

2 x iy5 

3x s y| 

4yi 

0 

2 x|y 5 

3x|y| 

4xsy| 5y5 



0 

-1 

0 

- 2 X 5 

-ys 

0 

7 2 

- 3 x 5 

-2x 5 y5 

-yi 

0 

- 4 x 5 

-3x|y 5 

- 2 x 5 y| 

-y| 

0 

-5X5 

-3x|y5 

- 2 x 5 y| 

-ys 

0 



1 

x 6 

ye 

4 

x eye 

v 2 

y 6 

V 3 

x 6 

x lye 

x ey| 

y| 

x 6 

x |ye 

2 2 
HYs 

x eyi 

ye 

V 5 

x 6 

x iyi 

x iye 

x eye 

yl 



0 

0 

i 

0 

x 6 

2y 6 

0 

4 

2 x 6 y 6 

3yi 

0 

x 6 

2x|y 6 

3x 6 y| 

4y| 

0 

2xgy 6 

3xgy| 

4x 6 y| 5y£ 



0 

-1 

0 

-2xg 

-ye 

0 

“3x| 

- 2 hy & 

-y| 

0 

-4x| 

-3x|y 6 

-2x 6 y| 

-ye 

0 

-5xg 

-3x|yg 

-2x 6 yg 

-ye 

0 



Bi(l,16) * -120 (Afj + A 13 A 21 - 0.5A u x 2 ) 


(25 j) 


B i(l,17) - ~12(A 11 A 3 2 + A 13 A 3£f + A 38 A 31 + A 39 A 33 + A^A^g 

(25k) 

+ A 15 A 21 - 0.5A 32 x 2 - A 31 xy - O.SA n y 2 ) 

B 1 (l,18) = -12(A 11 A 15 + A 13 A 25 + A 38 A 32 + A 39 A 34 + A lg A 31 

(251) 

+ a 15 a 33 - 0.5A 15 X 2 - A 32 xy - 0.5A 31 y 2 ) 

Bi(l,19) = -24(Aj5A 38 + A 25 A 3 g + A 16 A 32 + A 15 A 34 ~ A 15 x 7 

(25m) 

- 0.5A 32 y 2 ) 

Bj (1,20) = -120(A 15 A 16 + A 1 5 A 25 - 0.5A 15 y 2 ) (2Sn) 

Bl(2»7) = 6A 2 i (25o) 

B ].(2, 8) = 2A 33 (25p) 

B i (2,9) = 2A 34 (25q) 

Bi(2,10) = 6A 25 (25r) 

B l (2,11) = 24A 21 x (25s) 

Bj, (2jl2) = 6(A 33 x + A 21 y) ' (25t) 

B i (2, 13) = 4 (A 34 x + A 33 y) (25u) 

Bx (2,14) = 6 (A 25 x + A 34 y) (25v) 

B i (2 , 15) = 24A 25 y (25w) 

B 1 (2,16) = -120(Ai 1 A 2 i + A 23 A 2 i - 0.5A 2 iX 2 ) (25x) 
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Bi (2,17) = -12 (A 21 A 3 2 + ^■23^34 + A. 40 A 31 + Ajj.jA .33 + A 2 gAn 


(25y) 


+ A25A21 - O.5A34X 2 - A 33 xy - 0 . 5 A 2 iy 2 ) 

Bi( 2 , 18 ) = -I2CA21A15 + A23A25 + A4.0A32 + A41A34 + A 26 A 3 
+ A25A33 - O.5A25X 2 - A 3 4xy - 0 . 5 A 33 y 2 ) 

Bi ( 2 , 19 ) = -24 (A15A40 + A25A4.J + A26A32 + A25A34 - A 25 xy 

- 0.5A 3 4y 2 ) 

V 

Bi ( 2 , 20 ) = - 120 CA 15 A 26 + a| 5 - 0 . 5 A 25 y 2 ) 

\ 

An = -( J 11 D 11 + Ji2°13) 

A 12 = "( J 11 D 12 + J 12°23) 

A13 = -CJ11D13 + Jl 2 D 33 ) 

Ai 4 = -(J11D13 + J12D12) 

A 15 = _ C J 11 D 23 + J 12°22) 

A 16 = “( J 11 D 33 + <Jl2 D 23) 

A 2 I = "(^12^11 + ^22^13) 

A 22 = ”CJl2 D 13 + J 22 D 23) 

A 23 = “C J 12 D 13 + J 22 ^ 3 3 ) 

A 24 = "CJl2 D l3 + J22 D 12) 


(25z) 


(25aa) 

(25bb) 


(25cc) 


(continued) 
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A 25 = ~C J 12°23 + j 22 d 22) 


A-26 = "C j 12 D 33 * J 22 D 23) 

\ 

a 31 = A 14 + 2A 13 

a 32 = a 12 + 2A 16 

a 33 = A 24 + 2a 23 

a 34 = a 22 + 2A 26 

a 35 ~ a 33 + A 11 

a 36 = A 34 + A 31 ^ 

a 37 “ A 2S + A 32 

a 38 = a 13 + a 14 

a 39 = A 12 + A 16 

A40 = A 2 3 + A24 

A 4 I = a 22 + A 26 1 


(25cc) 

(concluded) 


7 . The matrix [b 3 ] relating (x 2 ) > the contribution of transverse shear 
to the vector of curvatures, to the generalized coordinates (m the equation 
fX 2 > = [B 3 ] {a}) is given by 


B 3 (l,ll) = -24A n 

(26a) 

B 3 (l,12) = - 6 A 31 

(26b) 

B 3 (l,13) = -4A 32 

(26c) 
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B 3 (l,14) = - 6 A 15 

C26d) 

B 3 (1,16) = -120A n x 

(26e) 

B 3 C1,17) = -12(A 32 x + A 31 y) 

(26f) 

= -12(A 15 x + A 32 y) 

C26g) 

B 3 C1,19) = -24A 15 y 

(26h) 

B 3 C2,12) = - 6 A 21 

(26i) 

B 3 (2,13) = -4A 33 

(26j) 

B 3 C2,14) = - 6 A 34 

(26k) 

B 3 (2,15) = -24A 25 

(261) 

B 3 (2,17) = -12(A 33 x + A 21 y) 

(26m) 

B 3 C2,18) = -12(A 34 x + A 33 y) 

(26n) 


B 3 (2,19) ~ -24(A 2 5X + A 34 y) 

(26o) 

B 3 (2,20) = -120A 25 y 

(26p) 

B 3 (3,ll) = -24A 21 

(26q) 

B 3 (3,12) = - 6 (An + A 33 ) 

(26r) 

B 3 (3,13) = “4(A 3 ^ + A 34 ) 

(26s) 

B 3 (3,14) = - 6 (A 32 + A 25 ) 

(26t) 

B 3 (3,15) = -24A 15 

(26u) 
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B 3 (3,16) = -120A 21 x 


(26v) 


8 3 (3,17) = -12 ^34 + A 3 i)x + (A 33 + An)y3 (26w) 
83 ( 3 , 18 ) = -12 f(A 2 5 + A 32 )x + (A 34 + A 31 )yJ (26x) 
83 ( 3 , 19 ) = -24 [a 15 x + (A 32 + A 25 )y] (26y) 
B 3 (3,20) = -120A 15 y (26z) 


where An , Aj 2 , . . ♦ , A34 are as given m equations (2Scc) . 

8. For plates with transverse shear flexibility the modified [hJ 
matrix, M , is given by subtracting the matrix M for each of the six 
grid points from the respective rows of a and 8 of the grid points m 

the M matrix. 

9. For plates infinitely rigid m transverse shear, 

[H 1 ] == [h] (27) 

10. The two constraint equations involving the coefficemts a 16 , 
a 17 , an , aig , and a 20 of the qumtic polynomial for transverse 
displacement so as to insure cubic edge rotation on the sloping edges of the 
triangular element are now entered as the 19th and 20th rows of M . 
i.e., the 19th and 20th rows are* 

'0 00000000000000 Sb^c (3b 2 c 3 -2b 4 c) (2bc 4 -3b3 c 2 ) (c s ~4b 2 c3) -Sbc 4 ' 

.000000000000000 Sa 4 c (3a 2 c 3 -2a 4 c) (-2ac 4 +3a3 c 2 ) (c s -4a 2 c3) 5ac 4 

This is now added as the 19th and 20th row of the £h'~| matrix to form 
[H'O matrix. [«"] is a 20 x 20 square matrix that is nonsmgular 1 

1 A numerical experiment to verify that H is nonsmgular for all practical 

element sizes is described m "New Triangular and Quadrilateral Plate-bending 

Finite Elements" by R. Narayanaswami , NASA TN D-7407, Apr 74. 
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11. The [h'O matrix is inverted. The first 18 columns of the 
[H"]- 1 matrix is denoted by the matrix £h , m ] (Size 20 x 18), i.e., « 

= The first 18 columns of jH'j] -1 (28) 

4.87.22.3 Stiffness Matrix Calculation for TRPLT1 (Subroutines KTRPLS and KTRPLD) 


The polynomial expressions, for variation of w and t within the 
element, are 


20 m. n 

L X 1 

ax y 
i=l 


(29) 


t 



r. s 

i i 

c.,* y 


(30) 


The values of m. , n , p. and q. are 
i * i r i n i 

mj = 0 ; m 2 - 1 , m 3 = 0 ; m^ = 2 ; m 5 = 1 ; m 6 = 0 , 

m 7 = 3 ; m 8 = 2 ; m 9 = 1 ; m 10 = 0 ; m n = 4 ; 

m 12 = 3 ; m 13 = 2 , mjij = 1 ; m 15 = 0 , m 16 = 5 , 

m i 7 - 3 , m 18 = 2 ; m 19 = 1 ; m 2 o = 0 

n i - 0 , n 2 = 0 ; n 3 = 1 , n 4 = 0 , n 5 = 1 , n 6 = 2 ; 

n 7 = 0 ; n 8 = 1 ; n 9 = 2 , n 10 = 3 ; n u = 0 , 

n i2 ~ 1 > n i3 ~ 2 > n i4 ~ 3 ; nj 5 = 4 , n 16 = 0 ; 

n 17 = 2 J n 18 = 3 , n 19 = 4 ; n 20 = 5 


(31) 


(32) 


r i - 0 ; 

r 2 = 1 , 

H 

to 

II 

o 

(33) 

</> 

II 

o 

V* « 

S2 = o ; 

s 3 = 1 

(34) 
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The coefficients a 3 to ago are generalized coordinates of the 
element and can be evaluated once the displacement vector is known. 

The coefficients cj , eg } and C3 can be evaluated from the 

specified thicknesses t* , t3 , and tj of the 3 corner grid points and 

the geometric dimensions of the element 


tja + t3b 
Cl = (a + b) 


(35) 


t 3 - t x 
° 2 ” (a + b) 


(36) 


c 3 = | (t 5 - cj) 


(37) 


where , t 3 , ts are evaluated from the values Ij , I3 and I5 
respectively. 

The elements of the symmetric portion of the stress-strain matrix 
eg are denoted by , ^12 > ^13 » Ggg > Gg 3 > and G 33 . 

A formula for the integral of the type x m y n taken over the area of 
the element is 


rr m n , . „ , v n+lr m+1 , , .m+l, m!n! 

JJ x y dx dy = F (m,n) = c {a - (-b) } ( m + n » 2 ) ; ( 36 ) 

The equation used is the stiffness matrix generation m generalized 
coordinates for plates infinitely rigid in transverse shear is given by 


( Vgen = 12 k £ S 1 %%% 


[g,,! m fm - 1) (m - 1) F(m + m + r, + r. + r. - 4, 

<- 11 1 j v 1 '"'•j •> v 1 j Icj k 2 k 3 

n + n + s. + s. + s-. ) + G 22 n n (n - 1) Cn - 1) F(m + m 
1 3 ki k 2 k 3 -' zz 1 1 J ^ j J v 1 j 


(39) 


+ r, + r, + r, , n + n + s. + s, + s. - 4) 
k l k 2 k 3 1 J k l k 2 k 3 


(continued) 
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+ (4G33® m n.n + Gi 2 {m n (m. - l)(n. - 1) + m.n. (m - 1) 
v 00 1 j 1 j i i 3 ] r j 


• (n - 1) ) F(m. + m + r, + r, + r, - 2, n. + n + s, 

^ i i 3 ki k 2 k 3 * i j ki 

+ s. + s, - 2) + 2Gi o{m.ra.n. (m. - 1) + m.n.m.(m. -1)} F(m 
k 2 k 3 ' l ] l i l j v j Ji ^ i 


(39) 

(concluded) 


+ m + r t + x, + r, - 3, n. + n. + s, + s, + s, - 1) 
3 kj k 2 k 3 n - v ^ v - 


“j % ' -k 2 ' -k 3 


+ 2G,-{m n.n (n. - 1) + m.n.n (n - 1)} F(m + m + r. + r. 
231 J l ] i u J 3 1 3 k x k 2 

* % - J > “i + n j * \ * % * % - 3) ] 


For plates with transverse shear flexibility, the expression for 
generalized stiffness matrix consists not only of the closed form expression 
but four additional integrals, given below, that are evaluated using numerical 
integration, i.e.. 


[K ]= & ](eq. 39) +//[ B 2 ] T [D] [b 3 ] dx dy 

6 s closed form 


*ff Cb 3 ] T [d] [b 2 ] dx dy + //&3] T 0>] [B,] dx dy (40) 
*//[» l] T 0= s ] [B,] dx dy 


0 >] 


matrix is obtained from the stress -strain matrix 



as 



c c.c, 

l J k 


W T k 

X J 



+s .+s 
J 


k ) 


(41) 


[G s ] = G t* 


1 

0 


0 

1. 


(42) 


where t* is the effective thickness for shear at the integration point and 
is evaluated from the user specified values TS1, TS3, and TS5. The 
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numerical integration formulae used are the seven-point integration 
scheme listed m Zienkiewicz 1 and are given below. 

For a triangle, the integrals of the form 

1 ,1-L 7 

I = / J fCLiL 2 L 3 D dL x dl 2 = £ w f. (LjLaLg) 

0 0 k=l 

where the points CLi, L 2 , L 3 ) and the weighting factors are as follows 


Point 

Triangular Coordinates 
l 3 

Weight , 2W^ 

1 

1/3, 

1/3, 

1/3 

0.225 

2 

«1 

6i 

Si 


3 

$1 

“i 

Si 

0.13239415 

4 

Si 

si 

«1 


5 

&2 

32 

s 2 


6 

S 2 


s 2 

0.12593928 

7 

32 

32 

a 2 




with 

otj = 0.05971588 gj = 0.47014206 

a 2 = 0.79742699 g 2 = 0.101286505 

The stiffness matrix m global coordinates is 

H»[E][T] T [H"0 T [ k ] gen [H...] [t] [e] T (43) 


1 O. C. Zienkiewicz, "Finite Element Method on Engineering Science," 
New York, London McGraw-Hill, 1971. 


133 



4. 87. 22. 4 Mass Matrix Calculation for TRPLT1 (Calculated m Stiffness 
Subroutines KTRPLS and KTRPLD) 

Two different mass matrices are used: the lumped mass and the 

consistent mass. The lumped mass matrix is calculated in the same manner 
as for TRIM6: 


m = -g- (pV + Ay) 


(44) 


where 


V, the volume of the element = cjF(0,0) + C2F(1,0) 
* C3F(0,1) 


(45) 


For each point, the diagonal mass matrix m element coordinates at 
all the grid points is 



0 

0 

m 

0 




i = 1, 2, 


6 (46) 


so that [M ] the element mass matrix has Cm. J matrices arranged 

O l 

diagonally. 

The mass matrix m the global coordinate system is obtained as 


[« gg ]=H [t] t [« ee ] [t] [e] T 


(47) 


If the parameter COUPMASS is set by the user, the consistent mass matrix 
will be formed. The jth element of the ith row of generalized mass matrix 
is given by 
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„ „ 3 m +m +r. n +n.+s, 

[mJ = p//E 3 ‘r 1 3 k tody 

J gen 


k=l 
m.+ra. n +n 


+ v 


ff x 1 ■* y, 1 ■* dx dy 


( 48 : 


= P E fc k F(m 4 * Hj + V ■>! * "3 - s k 5 ) 


+ y F(m + m , n •*• n ) 

1 3 1 r 


The mass matrix in global coordinates is 


(49) 


[m] .[b] [t] T [h- ] T [M ge J [h™] [t] [e] 


(50) 


4.87.22.5 Structural Damping Matrices for the TRPLT1 Element 
The structural damping matrices are 

&„] - S e bf,] (51) 

where g is the structural damping coefficient for the bending material 

C 

referenced. 

4.87 22,6 Stress and Element Force Calculations for the TRPLT1 Element 
(Subroutines STRP11 and STRP12 of Modules SDR2) 


1. STRP11 xs used to calculate the phase 1 stress-displacement 
relations . 

Frequent reference will be made to the equations from sections 4.87 22.2^ 
and 4 87.22.3 

The following data are calculated 

1 . [«■■■] - 20 x 18 Matrix relating generalized coordinates to grid 

point displacements. 
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2- &a] - 3 x 20 Matrix relating bending curvatures and generalized 

coordinates . 

3. [jB 3 - 2 x 20 Matrix relating curvature contribution of transverse 
shear strains and generalized coordinates 

4. [E] - element to basic coordinate transformation. 

5. [bj - 3 x 3 Matrix of elastic coefficients relating bending moments 
and curvatures. 

6. [G] -2x2 Matrix relating transverse shear forces and shear 
strains . 

- i = 1, 2, . . . ,6- Global to basic transformations. 

The following calculations are performed. 

[<] * H M [»"'] <‘- 2 > 

is a 3 x 18 matrix; this is 
as follows: 

r r * 1 * ! 

C S m 1 S j_ S M 1 j S M 2 j ’ * 

Each of the six matrix partitions is multiplied as follows : 

M = [\]M T W i - 1, 2, .... 6 (54) 

[ S D= BUMM C55) 

[Sg 3 1S a 2 x 18 matrix, this is split into six 2x3 matrix partitions as 
follows: 


split into six 3x3 matrix partitions 


I * *1 
! \\ 


(53) 





(56) 
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Each of the six matrix partitions is multiplied as follows: 


[e] T [tJ i = 1, 2, . . . , 6 (57) 

The 5x6 matrix [s^] is obtained as 

og= 

2. Phase 2 
(a) The vector of forces is computed as 

6 x 

£ [ Si ] (S9) 

where is the thermal moment vector If the thermal gradient is 

specified. 


/ M x V 

M 


M 


xy 


V 


M 


i “ 1 j 2 ^ ^ 6 


(58) 



(M t i ^ = -ogw; ( 60 ) 

where I is the moment of inertia of the cross section and T f is the 
i 1 

thermal gradient at vertex l of the element. 

The stresses and forces are evaluated at the vertices of the element, 
m addition, the stresses are also evaluated at the centroid. The simplification 
is made that the thermal moment vector at the centroid is the average of that 
at the vertices . 
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(b) With no given temperatures at the stress points, the stresses are 
then calculated from the equations 


M z 



\ 


M z 
y 

a y = I 


(61) 


M z 

- *z 


xy 


M x > and I for the 

appropriate points (vertices or centroid ) ; 
z values for the corner grid points are 
as those given in the PTRPLT1 card and for 
the centroid, the z values are the top 
and bottom fibre distances. 


If mean temperature 
three vertices. 


Tq and the gradient T' are specified at the 


\ 


" 



o 

X. 

1 

z 


M 

X 


K 

• - -T 

< 

M y 

> 

Q 

F 

N. 



M 

xy) 



+ &] {ot}T' - ~ (T q - T) [d] {a} i 


1 , 2 


(62) 


where T is the average temperature of the element. 

The principal stresses and angles are calculated using the same formula 
as for the membrane element TRIM6 (section 4.87.21 6). 

4.87.22.7 Thermal Load Calculations for the Bending Elements (Subroutine 
TL0DT1 , TL0DT2 and TL0DT3 of Module SSG1) 


The variation over the surface of the element of the mean temperature, 
Tq , and the thermal gradient at a cross section, T* , is assumed as a 
bilinear polynomial 

3 p q 

T 0 = E d * X Y 1 (63) 

i=l 


138 



(64) 


3 , P q 

T’ = E dj x v 1 

i=l 


so that the temperature at any point (x, y, z) is T = + T’ . 

The constants and d^ are evaluated from the values at the 

vertices and the reference temperature of the element 


dj = 


Tg^a + To3b 
(a + b) 


3 d 2 ~ 


t 03 - T 0 i 
"(a - + b) 


> d 3 ~ 7 - [ t 05 " di] 


( 66 ) 


dl - 


Tja + Tab 
(a + b) 


T3 - t} 


>^ = TTTbT J d 3 = F [Ts " d}] 


(67) 


It is convenient to define the elements of Tg 1 {o } as 

*- e J e 


G 11 = G ll« + G 12 a p + G 13 a p 
e l ©2 e 3 


( 68 ) 


G 22 “ G 12 a ei + G 22« e2 + G 23 a g 3 


(69) 


G33 = G 13 a ei + G 23 a e3 + G 33 a e3 


(70) 


The thermal load vector m generalized coordinates, (Fcren^ » will be 

evaluated m two stages, viz., the closed form expression {Fg en }l > due to 

the vector of curvatures in the absence of transverse shear and the 

numerically integrated expression {P t } 2 due to the contribution of 

gen A 

transverse shear to the vector of curvatures. The ith element of {P t }, 

gen 1 

is given by 


■ A 1, i, 1 1, v.,v! 


x (D^ - 1) + r x 


+ r 


10 + r i, + P, ' 2 < n , + s , + s + s, + q ) + G2 2 n.(n - 1 ) F(m 

x 2 x 3 3 1 x z x 3 3 1 1 1 


(71) 


+ r + r + r 
M 12 13 


+ p , n + s. + s. + s + q - 2) 
J 1 *1 12 13 J 


(continued) 
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+ G33 m-n, F(ra + r * r. + r. + p. - 1 , n + s + s 

i i i ij k 2 k 3 r j , i ij x 2 


* S 1 3 * q 3 ' 1)] 


( 71 ) 

(concluded) 


The load vector ( p g en ^2 is evaluated using numerical integration of 
the following expression. 


{pj en } 2 - 12 ff C G J { “e } T ' t3 dx iy t72) 

The generalized thermal load vector is 

(P t } = {P* >i + i? t } 2 (73) 

gen gen 1 gen * 

The thermal load vector m global coordinates is 

fp'lg - DO [t] T [h"'] T {P* en > C74) 

The forces are placed in the PG load vector data block. 
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4.87.23. TRSHL. Shallots Shell Triangular Element 


4.87.23.1 Input Data for TRSHL Element 


1. EST entries for TRSHL are 
Symbol 


EID 

SIL1, . . . SIL6 
0 

Mat ID 

m 


Ti, T 3 , T 5 

Mat ID, 

D 

1 1 > I 3> I 5 
Mat ID 

s 

TS1, TS3, TS5 
y 

Zll, Z21, Z13, Z23 
Z15, Z2S 
N I 

l 

X 

1 ' 1 = 1 ,.. 

Y 

l 

Z 

i 

TE1, TE2, . . TE6 


Description 

Element Identification Number 
Scalar indices of connected grid points 
Anisotropic material orientation angle 
Material-Identification Number for membrane 
behavior 

Membrane thickness at corner grid points 
Material-Identification Number for bending 
Area moments of inertia at comer grid points 
Material -Identification Number for transverse 
shear 

Thickness for transverse shear at corner 
grid points 

Nonstructural mass per unit area 
Distances Z1 and Z2 for stress 
calculations at three corner grid points 
Local coordinate system numbers 
and location of coordinates m the 
basic system for the connected grid 
points 

Element temperature at the six grid points 
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2. Coordinate system data 

The numbers , and Z ^ are used to calculate the 3 by 3 

basic-to-global coordinate transformation matrices T for points 
1 = 1, 2, 3, 4, 5, 6, (via subroutines TRANSD or TRANSS)-. 


3. Material data 
Symbol 

16 ] 


For mat. 
ID 

m 


P 


a 


x> 



xy 



V 


o 

s 


For mat. 


D 


Description 

3x3 stress-strain matrix 
Mass density 

Thermal expansion coefficients 
Reference temperature 
Structural damping coefficient 

Stress limits for tension, compression and shear 
3x3 bending stress-strain matrix 


ID b 1 

For mat. ( G Shear coefficient 

I s 

ID I 

s ' 

4.87.23.2 Basic Equation for TRSHL 


The calculations for the TRSHL element are very similar to those of TRIM6 
and TRPLT1 (sections 4.87.21.2 and 4.87.22.2 respectively) that only the essential 
details are given here. 


The displacement transformation matrix from basic coordinates to m-plane 
coordinates is 



il 

Jl 

k l 

0 

0 


3-2 

12 

k 2 

0 

0 


13 

13 

k 3 

0 

0 


0 

0 

0 

il 

11 


0 

0 

0 

12 

5 2 


0 

0 

0 

13 

J3 


( 1 ) 
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where 1 , j and k are as given by equations (4) , (5) and (6) of section 
4.87.22.2. No transverse shear effects are considered for the TRSHL element. 

The matrix [H] relating grid point displacements and the generalized 

coordinates (m the equation {vj = [H] {a} ) is similar to that for TRIM6 

and TRPLT1. The inverse matrix relating generalized coordinates to the grid- 
point displacement vector is a 32 x 30 matrix and is given by 


M - 1 


Hf 1 

1 

| 

1 

0 1 

6x6 




1 

| 

Hf 1 1 


1 

6x6 J 

- 

— 1 - 



I 

i H m 


1 

i 

* 20x18 


( 2 ) 


where Hi is a 6 * 6 matrix given m equation (21), section 4 87.21.2 and 
H ,M is a 20 x 18 matrix given in equation (28) of section 4.87 22 2. 


4.87.23.3 Stiffness Matrix Calculation for TRSHL (Subroutine KTSHLS and KTSHLD) 


The polynomial expressions for variation of u , v , w and thickness t 


within the element are 





32 


m. 

n 


U - 

£ 

a 

l 

X 

l 

Y 

(3) 


1=1 





32 


P, 



V - 

£ 

b 

i 

X 

y 1 

(4) 


1=1 






32 


V 

S 


w - 

£ 

C 

1 

X 

y 1 

(5) 


1=1 
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t = 
m 

3 

Z 

i=i 1 

t. u 
x 1 y 1 


(6) 

\ - 

s 

E ^ 

l = l 

t* u! 
X 1 y 1 



The values of nu , 
are 

n i ’ Pi > % • V 1 ’ s i 1 

h * U i * 

t ' and 

i 

V 1 

= 1,32): 

0, 1, 0, 2, 1, 0, 26*0 



V 1 

= 1,32): 

0, 0, 1, 0, 1, 2, 26*0 



P x Ci 

= 1,32): 

6*0, 0, 1, 0, 2, 1, 0, 20*0 



q t Ci 

= 1,32): 

6*0, 0, 0, 1, 0, 1, 2, 20*0 



v ± Ct 

= 1,32): 

12*0, 0, 1, 0, 2, 1, 0, 3, 2, 

1, 0, 4, 

(7) 


V 

3, 2, 1, 0, 5, 3, 2, 1, 0 



s id 

= 1,32): 

12*0, 0, 0, 1, 0, 1, 2, 0, 1, 

2, 3, 0, 



1, 2, 3, 4, 0, 2, 3, 4, 5 

t.(i = 1,3): 0, 1, 0; t^:0, 1, 0, s^O, 0, 1; s^O, 0,1 

The coefficients a^ , fa^ and c^ are undetermined parameters such that 
a x = 0 i = 7 to 32 

b> i = 0 l = 1 to 6 and 13 to 32 (8) 

c = 0 i = 1 to 12 

l 

a^, l = 1 to 6; b^, i = 7 to 12, and c^, l = 13 to 32 can be determined once 
the element displacement vector is known. 
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The coefficients d and d' , i = 1 to 3 are evaluated from the 

11 

user specified values of the membrane thickness and the area moments of inertia, 
respectively, by equations similar to those for TRIM6 and TRPLT1 elements. 

The equation used m the stiffness matrix generation m generalized 
coordinates is (following the procedure outlined m sections 5.8.6 and 
5.8. 7), the jth column of the ith row of the generalized stiffness matrix ’is 
obtained as 


K. . = T ("Giifm m. d. F(m. 
ij L Ll \i j k ^ i 


- h., m 1 d k FCm L + + t k - 1, n 4 ♦ 

- hi* d k PfMj + r i + t k - 1, + s. + u k ] 

* h? d k F( Tl + * t k , » t + Sj - u fc )) 

* G 22 ( qi q. d k Ftp. + Pj + t k , q, * » u k - 2) 

- h 6 d k F( Pl * ♦ t k , ^ * s. - u k - 1) 

- h 6 d k F(r x ♦ P, ♦ V S i * + “k ' 1] 

. h| d k Ft^ . rj * t k , \ + u k )) 

+ G33 (n n d, F (in + m + t. , n + n. + u, - 2) 

^ \ i j k v i j k*i j k 

+ n. p. d. F(m + p + t, - 1, n + q + u, - 1) 
i k v i k * i k 

- h 5 d k F(m x + r^ + t k , tl % + s^ + u R - 1] 

. Pj "j d k n Pl * »j - \ - 1. q k + ♦ \ - 1) 

♦ Pi Pj d k p( Pi + Pj * \ - 2 - <k + i, + V 

- h 5 P 1 d k F CPl + r, * t k - 1, q x ♦ s 3 * u k ) 

- h 5 n d, F(r + m + t. , s + n + u, - 1) 

j k l ■) k l j k 


( 9 ) 


(continued) 
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- h s p 3 d R F(r. ♦ Pj ♦ \ - 1. »i ♦ <Ij ♦ \) 

* hi d k FC*i - r. * t k , 5i *y u k )J 

+ Gl2 ( m l \ F( \ Pj + C k - U n i + 1 j + u k - 1J 

- h 6 m i d k F(m i + + t k - 1, n. + s. + UjJ 

“ h 4 d k F(r^ + Pj + t k , s j , + q^ + u R - 1) 

+ 2h 4 h 6 d k + t k , s. + s.. + 

+ \ m j d k F ^i + m 3 + h~ '> \ + n j +u k~ « 

- h4 q ± d k F(p i + + t R , q ± + s. + ^ - 1) 

- hg nu d k F(r ± + in,. + tj. - 1, s A + 

+ Gi3^. n^ d k FO^ + + t k - 1, n x + + u k - 1) 

+ m x Pj d k F(m i + p 3 * h - 2 ’ n i + *3 + V 

- h 5 d k F (nn + r + - 1, n. + s. + 1^) 

- h 4 d fc F(r i + + t k , s 1 + + u R - 1) 

- h 4 Pj d k F(r. + P-, + \ - 1, q ? + u^) 

+ 2h 4 h 5 d k F(r 2 + + t^, s i + s. + u k ) 

+ n. m. d k F^ + m. + - 1, 1^ + n.. + u k - 1) 

- h 4 n. d, F(m + r + t, , n. + s + u, - 1) 

H 1 k v 1 j k’ 1 j k J 

+ p x m j d k FC Pl - m. * t k - 2, q. + + u R j 

- h 4 Pi d k F( Pi + r j + t k - 1, q x + + u k ) 

- h 5 Bj d k F( Tl ♦ m. ♦ t k - 1, Sj ♦ n 5 + u k )) 

♦ 023^ d k F(p. ♦ m + t k , + n 3 + u k - 2 ) 

+ q i Pj d k F(p i + P 3 + l k ' X > q i + <1] * u k - 15 


(continued) 


146 



• h5 % ‘‘k Ptp i * r J + V ’i + + u k - 15 

- he n d, F(r + m. + t, , s + n. + u, - 1) 

D ik^i j k’ 1 3 k J 

- h 6 Pj d k FC^ * Pj * t k - 1, s . * q . * u k ) 

+ 2h 5 h 6 d fc F (r^ + + t k , s i + + u k ) 

- n x q } d k F( Bi * Pj * V \ ♦ q., - \ - V 

- h s n. d k F( Bi * r., * t k , n a + s j + d k - 1) 

+ P i % d k F(p i * P 3 * l k - *• <>i * < 1 ; * »k - 15 

- h 6 Pl d k F(p. ♦ ♦ t* - 1. q 1 ♦ Sj + u k ) 

- hs q ? d k F(r x * p. * t k , s. * qj + u k - 1))] 


3 3 3 

+ Z L £ C T2 d L d L d L CG n r ! r q ^ - X ^ r q - u 

k 1= l k 2 =l k 3 =l 12 K 1 k2 k 3 1 J 1 J 

* F(r. + r + t ' + t' + t.' - 4, s + s + u' + u,’ + u,' ) 

i 3 ki, k 2 k 3 l 3 k x k 2 k 3 J 


C9) 


+ G 22 s. s (s - l)(s - 11 F(r + r + t ' + t ' + t ’ (concluded) 

i J i 3 l j ki k 2 k 3 


+ s. + 

l 


+ u, 


+ u: + u 


ki k 2 


ks 


4) 


+ (4G 33 t i t, + Gi 2 {r x (r. - 1) (s^ - 1) 

♦ r. Sj (r. - l)( Si - 1)}) Ff^ * r 3 + * t£ $ 


2 , 


+ s + * 
l 

5 + 

J 

u ki * 

U» 

k 2 


- 2) 




+ 2G 13 {r i 

r i 
J 

5 j cr x 

- 1) 

+ r. 
i 

r j S x tr j 

- I)) F(r x 

+ r 

3 

+ t,' + 

ki 

t 1 

k 2 

+ K 

k3 

- 3, 

s + 

l 

S J * ”ii * 

“i 2 


- 1) 

+ 2G 23 {r^ 

s * 

i 

3 (S 
3 1 

- 1) 

+ r 

l 

s d s j (S J 

- D| 

F(r i 

+ r. 

3 

+ t,' + 

*1 

t k 2 

+ t,' 
*3 

- 1, 

s. + 
i 

s + uJ + 
3 ki 

u’ 

k 2 

*“*3 

~ 3)] 
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The generalized stiffness matrix can be transformed to the element and global 
coordinates by transformations similar to those for TRIM6 and TRPLT1 elements. 

4.87.23.4 Mass Matrix Calculation for the TRSHL Element (Calculated in the 
stiffness _ subroutine KTSHLS and KTSHLD) 


Two different mass matrices are calculated: lumped mass and consistent 

) 

mass. The calculations are the same as for TRIM6 and TRPLT1. 

4.87.23.5 Structural Damping Matrices for the TRSHL Element 

The calculations are similar to those for TRIM6 and TRPLT1 elements. 

4.87.23.6 Stress and Element Force Calculations for TRSHL Element (Subroutines 
STRSL1, STRSLV and STRSL2 of module SDR2) 


The calculations are similar to those of TRIM6 and TRPLT1 elements. 

4.87.23.7 Thermal Load Calculations for the TRSHL Element (Subroutines TLgSDSL 
of module SSG1) 


The calculations are similar to those for TRIM6 and TRPLT1 elements . 

4.87.23.8 Differential Stiffness Matrix Calculations for the TRSHL Element 
(Subroutine DTSHLD of module SSG1) 


The steps leading to the calculations of the differential stiffness matrix 
are given m section 7.3.6 of the theoretical manual (pages 66 to 70 of this 
report) . 
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APPENDIX D 

Updates to the NASTRAN Demonstration Problem Manual 
for the addition of TRIM6, TRPLT1 and TRSHL elements 
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Demo Problem 1.3-4 (a) 


Analysis of a Free Rectangular Plate with thermal loading using higher 
order triangular membrane TRIM6 element The quarter section of the plate 
shown m figure 1, is discretized using TRIM6 element. Discretization is given 
on page 1.3-4(b), figure 3(a). 

The graphs for the measured stresses and at x = 1.5 shown on 

pages 1.3-5 and 1.3-6. The results obtained by this analysis are not included 
m the same graph, since for the chosen mesh the stresses are evaluated at 
locations different from those shown in the graph. However, good agreement 
is seen for the stresses for the chosen mesh. 


1 3-4 (a) (1/1/77) 




BLAK& NOE 
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Demo Problem 1.7-5 


Triangular Shallow Shell Element 

Two problems, (1) that of a spherical cap, and (2) that of a cylindrical 
shell roof, are considered. These are the same two example problems analyzed 
m reference 23. 

(1) The spherical cap, with finite element discretization is shown m 
figure 5. Due to symmetry, only one fourth (quarter] of the cap was analyzed. 

Good agreement m deflections at the center of the cap is obtained even 
with relatively coarse mesh sizes as shown m Table 1. Even though the results 
appear to be oscillating about the exact value, the percentage error in the 
converged solution is very negligible. 

(2] The shell is shown in figure 6 along with pertinent dimensions and 
associated material properties. The finite element discretization for the shell 
is shown m the same figure. Due to symmetry, only one fourth (quarter) of 

the shell was analyzed. 

Results for the shell roof problem and the exact solution reported by 
Cowper et al (ref. 23) are given m Table 2. Reasonable agreement is seen 
between the finite element and the exact solutions. 

The values given m Table 2 are obtained from a stand-alone program wherein 
the global stiffness matrix had 5 d.o.f. per grid-point, viz., u , v , w , 
a and 3 . This is consistent with shallow shell theory. In NASTRAN, the 
global stiffness matrix has 6 d.o.f. per grid point, viz., u , v , w , a , 

3 and y . It is necessary therefore to constrain the sixth degree of freedom 
at all grid points where all the elements connected to that grid point are in 
the same plane. This requirement is to ensure that the global stiffness matrix 
is nonsingular for a given sufficiently supported condition of the structure. 
Theoretically, however, the above requirement is equivalent to the introduction 
of additional constraints on the problem and hence the solution obtained from 
NASTRAN will be lower bounds to the actual values obtained from the stand alone 
program and given m table 1. The values obtained from NASTRAN and CTRSHL elements 
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are given m table 2. For shells that axe strictly shallow, the solution from 
NASTRAN will approach that obtainable from stand alone programs based on a strict 
application of shallow shell theory. 

An alternative to solve problems where shell is only marginally shallow, 
as the example discussed herein, is to use combination of TRIM6 and TRPLT1 
elements. The result of using a 2 x 4 and 3x3 mesh of CTRIM6 and CTRPLT1 elements 
from NASTRAN is given in table 4. Note that the values are very close to the 
exact values. 


1.7-6 (1/1/77) 
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Table 1. Center deflections for spherical cap problem. 


Etw 

Deflection 6 = — 

C PR 2 


Finite 

Element 

Grid 

— = 0.02 
L 2 

— = 0.005 
L 2 

1 x 1 

1.15107 

1.13951 

CM 

X 

CM 

1.00774 

0.99178 

3x3 

1.00452 

1.00177 

4 x 4 

1.00437 

1.00084 

Exact 

1.00978 

1 . 00043 


1.7-7 Cl/1/77) 
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1.7-8 (1/1/77) 


Cn 

ON 


Table 2. Results for a cylindrical shell roof from stand-alone program (values are m element 
coordinates) . 


Finite 

Element 

Grids 

10u. 

A 

(m.) 

W B 

(in.) 

10v B 

(in.) 

10 ”c 
(in. 3 

10" 3 N „ 

xxB 

(lb./m.) 

10“ 3 M _ 

yyC 

(lb. m./m.) 

10 ' 2M X xC 
(lb. in. /in.) 

1 x l 

-0.45168 

-0.29100 

-2.48424 

-4.0700 

2.4659 

, 0.7685 

2.8520 

2x2 

-0.7812 

-12516 

-4.77312 

-2.1344 

4.2801 

-0.9395 

-0.8896 

3x3 

-1.09590 

-2.49876 

-7.12872 

-1.3606 

5.4948 

-2.0283 

-1.1136 

4x4 

-1.2939 

-3.4332 

-8.57580 

2.2224 

6.0277 

-2.3828 

-1.7912 

2x4 

-0.9041 

-2.2815 

-6.15 

0.888 

5.1312 

-1.415 

-2.0196 

3x6 

-1.1244 

-3.6227 

-8.5968 

3.1031 

6.1862 

-1.9414 

-1.8912 

4x8 

-1.3845 

-4.1526 

-9.5295 

3.9238 

6.4839 

-2.0459 * 

-1.6724 

5x5 

-1.4160 

-3.88152 

-9.29000 

2.8182 

6.3279 

-2.3538 

-1.9770 

6x6 

-1.4733 

-4.09176 

-9.76992 

3.0900 

6.4444 

-2.3242 

-2.0638 

Exact 

-1.51325 

-4.09916 

-8.76147 

5.2494 

6.4124 

-2.0562 

-0.9272 




Table 3. Results for cylindrical shell roof from NASTRAN using 
CTRSHL elements (values are in global coordinates) . 


Finite 

Element 

Grids 

U A (in.) 

lV g (in.) 

V B (in.) 

\ (m.) 

2x4 

-0.0945 

-1.6437 

-0.5181 

0.1938 

3x3 

-0.09054 

-1.7309 

-0.4801 

0.3813 

Exact 

-0.151325 

-3.70331 

-1.96372 

0.52494 


1 7-9 (1/1/77) 
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Table 4. Results for cylindrical shell roof from NASTRAN using 
CTRIM6 and CTRPLT1 elements (values are m global 
coordinates) . 


Finite 

Element 

Grids 

U A (in.) 

W fi (in.) 

V B (in.) 

W c (in.) 

2x4 

-0.1233 

-3.4162 

-1.7445 

0.4287 

3x3 

-0.1335 

-4.2560 

-2.1226 

1.1007 

Exact 

-0.151325 

-3.70331 

-1.96372 

0.52494 
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DIAPHRAM 


DIAPHRAM 


FREE 

EDGE 


L = 50 ft. 
(15.24 m) 


FREE 

EDGE 



E = 3.0 x 10 6 psi (2.068 x 10 11 N/m 2 ) 
v = 0.0 
t = 3.0 xn. (7.62 cm) 

shell weight = 90 lb./sq. ft. (4.309 x 10 3 N/m 2 ) 


Figure 6 . Geometry of cylindrical shell roof and finite element 
idealization 


1 9-12 (1/1/77) 
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Demo Problem 1.8-5 


Analysis of a Beam Using TRIM6 


A. Description 

The cantilever beam shown on page 1.8-4 is modeled with the NASTRAN TRIM6 
element as shown in figure 2, page 1.8-7. This problem demonstrates the 
analysis of a beam subdivided into the six noded triangular membrane elements. 

The loads were chosen to approximate the stress distribution due to a 
moment on one end of a beam. The other end is constrained to resist the moment. 
The plane of symmetry is not used. 

B . Input 

1. Parameters Similar to those listed on page 1 8-1. 

2. Boundary Constraints, on x = 0 plane, u x = u^ = 0 . 

3. Loads- total moment = M^, = 2.048 x 10 3 . This moment will produce 

bending about the z-axis. It is modeled by a set of axial loads at 

x = Z which, m turn, represents an axial stress distribution: 

a = 1.5 . 

x y 

4. Subcase - 1 Consistent loading. 

5. Subcase - 2: Lumped loading. 

Considering strip near extreme fiber 

F = -J- (1. 5 x 6 + 1.5 x 8) 2 4 = 84 

x 2 

C. Analysis and Results 

Analysis, refer to pages 1 8-1 and 1.8-2 
Results 
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Comparisons of Displacement 


Grid Pts. 

Theory (10"^) 

Mx^ 

7 " 2El 

Consistent Loading 
Subcase 1 (10”^) 

Lumped Loading 
Subcase 2 (10 -tf ) 

3 

0 

0 

0 

13 

.0625 

.0515 

.0523 

23 

.25 

.2377 

.2467 

33 

.5625 

.549 

.5744 

43 

1.0 

.985 

1.0172 


Comparisons of Stress 


Figures 3(a) 
Referring to 
Grid point 5 
Grid point 3 

Grid point 1 


and 3(b) show stresses obtained from the analysis, 
figure 3a, subcase 1, we have stress at 
= 9.567 (C) 

1.64 + 1.31 + 1.31 . . 

3 1-4 (C) 

_ n _ 1 . 7 - * ._ 1 . 5 . 1 46 _ 13>S8 (T) 


Referring to 
Grid point 5 
Grid point 3 

Grid point 1 


figure 3b, subcase 2, we have stress at 


= 9.87 (C) 

1.64 + 1.35 + 1.35 


= 1.48 (C) 


12.1 + 15.96 


= 14.03 (T) 


If the cantilever beam is discretized with the same type of mesh and the same 
number of elements, hut the diagonal oriented m the opposite direction to 
figure 2, i.e., as shorn m figure 4, then the stresses for subcase 1, at 
grid points 5, 3, and 1 would be 

Grid point 5 = 13.58 (C) 

Grid point 3 = 1.4 (C) 

Grid point 1 = 9.567 (T) 

1.8-6 Cl/1/77) 
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Therefore the stresses at grid points 5, 3, and 1 are taken as the average 
of the two types of meshes (i.e., figure 2, and figure 4). 


Therefore, for subcase 1, we have stresses at 
9.567 + 13.58 


Grid point 5 = 


= 11.575 CCD 


Grid point 3 = — - * 1 - = 1.4 (C) 
Grid point 1 = 1 , 3-58 + 9.567 = 11>575 


Conclusion 


Stresses at Grid Point 

Theory 

NASTRAN 

TRIM6 

5 

12 

11.575 

3 

0 

1.4 

1 

12 

11.575 
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at y-8, a = 1 . 5y = 1.5 x 8 = 12 
X 

Total Axial Force = F^ = -|-x 8 xl 2 x 4 = 192 



x (1.5 x 2)4 = 12 


Considering strip m between N~A. nnd extreme 
fiber F = — (1 . 5 x 2 + 1.5 x 6)4 x 4 « 96 







Demo Problem 1.9-4 


Thermal and Applied Loads on TRIM6 Elements 

A. Description 

This problem demonstrates the use of the TRIM6 elements. Ten triangular 
membrane elements are used to model a 2 x 1 x 10 beam. The dimensions and 
boundary conditions are shown in figure 2. Two loading conditions are applied: 
axial stress and thermal expansion. Symmetry boundary conditions are used. 

B . Input 

1. Parameters, similar to those listed on page 1.9-1. 

2. Boundary Constraints: u^ = u^. = 0 at x = 0 , u^. = 0 at y = 0 . 

3. Loads: Subcase 1, consistent loading 

F = 24 x 10 3 (total axial force) 

x 24 

Total force on symmetric part = - 12 


1 x 12 
6 


1 x 12 
4 

Subcase 2, Thermal loading 

T = 60° (Uniform temperature field) 

T q = 10° (Refeaysa.ee temperature) 

6 ' 


4 

1 


Force divided into the ratio of 

, , _ 1 x 12 4 x 12 

1.4 1, i.e.j , } * , 


Subcase 3, Lumped loading 


1 

2 

— 1 


Force divided into the ratio of 

, . 1 x 12 2 x 12 

l z . 1 , i.e, , . , . » 
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C. Analysis and Results 

Analysis: refer to page 1 9-2 

Results : 



Subcase 2 


TRIM6 Sol. 


0.109 

0.2093 

0.3093 

0 4093 

0.5093 

0.6093 

0.7093 

0.8093 

0.9093 

1.00937 


Displacement (U ) 


Graph is given on pages 1.9-7 and 1 9-8 


Conclusion 


The results of all three subcases are exact to the single precision 


limits , 
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Demo Problem 1.11-4 (a) 


Analysis of a simply supported rectangular plate with a thermal gradient 
using higher order triangular bending TRPLT1 elements: The quarter section of 

the plate is discretized using TRPLT1 element. Discretization is given on 
page l.ll-4(b), figure 2. 

For input and theory, refer to pages 1.11-1, 1.11-2, and 1.11-3. 


Result : 


The maximum displacement obtained was 0.5935, a difference of about 5 
percent from the analytical value. A more refined mesh is likely to yield 
closer values to the exact ones. 


The graph for the moments M , M and M obtained by analysis 
r x ’ y xy 1 

at x = 0.5 is shown on page 1.11-6. The results obtained by this analysis 
is not included m the same graph, since for the chosen mesh the moments are 
evaluated at ^locations different from those shown m the graph. However, good 
agreement is seen for the moments for the chosen mesh. 
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Figure 2. Model of simply supported rectangular plate using 
TRPLT1 element. 
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Demo Problem 1.11-8 


Deflection of Thick Rectangular Plate Using TRPLT1 Element 

The TRPLT1 element is used for solving moderately thick to thick plate where 
the effects of transverse shear are present. 

The rectangular plate m figure 1, shown on page 1 11-11, is discretized 
by using higher order triangular bending element TRPLT1. Because of symmetry, 
the quarter section of the plate is discretized and details of the discretization 
are given m the same figure. 

Four different mesh sizes are used for Q-mesh and P-mesh The result is 
tabulated for the simply supported and clamped edge with two different ~ ratios, 
on page 1.11-9 and 1.11-10 respectively. 


Input : 

E = 3.0 x 10 7 lbs/m. 2 
v = 0.3 

q = 1000 lbs/m. 2 

q = 1000 lbs. 

* = 2 
a 

t = 1.0 m. 


(Young's modulus) 

(Poisson's ratio) 

(Uniform distributed load) 
(Concentrated load at center) 
(Length/width) 

(Thickness of the plate) 
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Table 1. Central deflection of simply supported rectangular plates 
(Q-mesh) including the effects of transverse shear. 


Number 

of 

Elements 

per 

Side 

N 

Concentrated Load 
at Center 

Uniformly Distributed 
Load 

t = 100 

t-« 

§- = 100 
t 

1=4 

2 

21.3373 

22.22 

10.983 

11.2651 

4 

17.7854 

20.0133 

10.2084 

10.5612 

8 

16.9276 

19.4859 

10.1396 

10.4831 

12 

16.73 

19.5322 

10.1330 

10.5468 
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Table 2. Central deflection of clamped rectangular plates (Q-mesh) 
including the effects of transverse shear. 


Number 

of 

Concentrated Load 
at Center 

Uniformly Distributed 
Load 

-Ci JL wlllLsXl 

per 

Side 

N 

f - 100 

t = 4 

T 35 100 

f-4 

2 

10.4330 

12.5582 

3.942 

4.5086 

4 

8.4230 

10.565 

2.7932 

3.133 

8 

7.6283 

10 0845 

2.5953 

2.9375 

12 

7.4336 

10.1791 

2.561 

2 9475 
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Figure 1. Discretization and schedule of rectangular plate. 
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Demo Problem 1.11-12 


Analysis of Rectangular Plate Using TRPLT1 Element 

This problem demonstrates the accuracy of TRPLT1 element m the evaluation 
of moments in rectangular plate. 

The rectangular plate with discretization, shown on page 1 11-9 (figure 1) , 
is analyzed. 

Two different mesh arrangements are used, i.e., Q-mesh and P-mesh. The 
result is tabulated for the points marked on x~ and y-axis, as shown in 
figure 2, page 1.11-17. 

Input : 

E = 3.0 x io 7 lbs/m. 2 
v = 0.3 

q = 1000 lbs/m. 2 

5 -= 2 

a 

t = 1.0 in. 

N = 12 
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Table 1. Numerical factor 8’ m the equation M x = 8’ q a for bending moments 
of simply supported rectangular plate under uniformly distributed load 
(along y = 0) . 


X 

NASTRAN - TRP LT 1 

Exact 

8* 


8’ (y = 0) 


Q-mesh 

P-mesh 

Average 

0.5a (center) 

0.1001332 

0.1022486 

0.1012 

0.1017 

0.416667a 

0.105326 

0.1070887 

0.10620733 

0.098875 

0. 33333a 

0.099048 

0.10247159 

0.100759796 

0.09075 

0.25a 

0.087942 

0.0923552 

0.0901486 

0.076875 

0.166667a 

0.0710574 

0.0751297 

0.07309355 

0.058 

0.083333a 

0.04853 

0.0393438 

0.0439369 

0.032 

0a (end) 

0.022253 

0.0254971 

0.02387506 

0.0 
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Table 2. Numerical factor £>{ m the equation My = q a for bending moments 
of simply supported rectangular plate under uniformly distributed load 
(along y - 0) . 


X 

NASTRAN-TRPLT1 

Exact 

Pi 

(3* (y = 0) 

Q-mesh 

P-mesh 

Average 

0.5a (center) 

0.046605 

0.0470838 

0.0468444 

0.0464 

0.416667a 

0.047194 

0.0474771 

0 0473355 

0.045 

0.33333a 

0.0432597 

0.0441237 

0.04369171 

0.0410625 

0.25a 

0.0369394 

0.0379362 

0.037438 

0.0344375 

0.166667a 

0.0280123 

0.02855153 

0.028282 

0.0252 

0.08333a 

0.01630765 

0.0100907 

0.0132 

0.013875 

0a (end) 

0.006676 

0.01902541 

0.0128507 

0.0 
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Table 3. Numerical factor $" m the equation M x = g" q a for bending moments 
of simply supported rectangular plate under uniformly distributed load 
(along x = 0) . 


Y 

NASTRAN-TRPLT1 

Exact 

S" 


6" (x = 0) 


Q-mesh 

P-mesh 

Average 

0.5b (center) 

0.1001332 

0.1022486 

0.1012 

0.1017 

0.416667b 

0.103584 

0.1087446 

0.1061643 

0.099625 

0.33333b 

0.0969304567 

0.1131157037 

0.1050231 

0.09025 

0.25b 

0.0850907 

0.1218893745 

0.10349 

0.080625 

0.1666667b 

0.0674407 

0.14173259 

0.1045866 

0.06175 

0.08333b 

0.04365204 

0.1811511029 

0.1124016 

0.033625 

0 . 0b (end) 

0.01275513 

0.037393 

0.025074 

0.0 
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Table 4. Numerical factor By m the equation My = 6" q a for bending moments 
of simply supported rectangular plate under uniformly distributed load 
(along x = 0) . 


y 

NASTRAN-TRPLT1 

Exact 

B’i 

ey (x * 0) 

Q-mesh 

P-mesh 

Average 

0 . 5b (center) 

0.046605 

0.0470838 

0.0468444 

0.0464 

0.416667b 

0.0486921 

0.05105566 

0.049874 

0.0465625 

0.33333b 

0.049534 

0.057433 

0.0534835 

0.0465625 

0.25b 

0.049873 

0.0679783 

0 05892565 

0.0460625 

0.166667b 

0.04794 

0.0850651 

0.06650255 

0.041125 

0.083333b 

0.03933845 

0.109444 

0.0743913 

0.0288125 

0.0b (end) 

0.001104162 

0 01272338 

0.0069137754 

0.0 


1.11-16 (7/1/76) 








Demo Problem 3.1-6 


Vibration of Tapered Rectangular Plates 


This problem demonstrates the use of the higher order triangular bending 
element CTRPLT1 to solve problems in vibration of thin isotropic plates. 

The structural problem consists of a linearly tapered rectangular plate 
with two different support conditions, namely ( 1 ) simply supported and 
(ii) cantilever. 

( 1 ) Linearly tapered simply supported rectangular plate The model as 
shorn in figure 4a uses only half of the plate due to symmetry. The plate 
thickness is given by: 


t 



(1 


+ K - ) 

a 


CD 


where k is a constant determining the rate of taper. Two different mesh sizes 
of the finite element model, 1x2 and 2x4, are used. Nondimensional 
fundamental frequencies for rectangular plates for three different aspect 
ratios -g- and k = 0.5 and 0.8 are presented m table 2. 

The frequency parameter is defined as* 


ft = 



C2) 


where to is the circular frequency, a is the length, p is the mass density, 
t Q is thickness and D q is the bending rigidity. Analytical results from 
reference 20 are also shown for comparison. 

(li) Linearly tapered cantilever rectangular plate The plate is 
idealized with a mesh size of 2 x 4 or 16 elements, as shown m figure 4b. 
Results of frequency parameters ft^ as defined m equation (2) , where m 
and n represent the number of nodal lines perpendicular and parallel to 
the support, respectively, using TRIA2 and TRPLT1 are shown m table 3 
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Constant thicknesses of 0.0405 in. CO. 1029 cm) and 0.1215 in. (0.0386 cm) were 

/ * 

used when modeling with TRIA2 element. Experimental data obtained by Plunkett 
in reference 21 are also given. 

Tables 2 and 3 show that very good results have been obtained using the 
higher order plate element. 
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Table 2. Fundamental frequency for linearly tapered rectangular 
plates simply supported on all edges* v = 0.3. 


Aspect 

Ratio 



NASTRAN 

TRPLT1 

Finite 

Element 

Layout 


1x2 

Theory 


1x2 

2x4 

Theory 


1x2 

2x4 

Theory 


Frequency Parameter 
/Pt n \ 1/2 

a ■ (-of) 

Taper Rate 

Taper Rate 

k - 0.5 

k = 0.8 

14.662 

16.242 

15.304 

16.994 

24.171 

26.901 

24.454 

— 

24.556 

27 354 

58.560 

60.346 

64.770 

60.982 

67.500 
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Table 3. Frequency parameters for a linearly tapered rectangular 
cantilever plate; v = 0.3. 


Mode 

Frequency Parameter = 

mn 

2 / pt 0\ 1/2 
“«m a ( D 0 j 

NASTRAN 

Experiment 

m 

n 

TRIA2 

TRPLT1 

0 

0 

2.28 

2.25 

2.47 

1 

0 

9.8 

10.0 

10.6 

0 

1 

14.5 

13.6 

14.5 

1 

1 

23.8 

27.0 

28.7 

0 

2 

/ , 

35.9 

32.8 

34.4 

0 

3 

51.5 

47.3 

47.4 

2 

0 

31.0 

53.3 

52.5 

1 

2 

64.0 

57.7 

54.0 
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Demo Problem 5.1-5 


Buckling of Columns and Plates 


The out-of-plane buckling of plate elements is evaluated from the 
differential stiffness matrix of bending plate element TRPLT1 due to membrane 
prestress effects obtained from a membrane analysis using TRIM6 elements. 

To solve out-of-plane buckling of plates, a membrane-bending combination 
element is necessary. TRSHL is such a combination element with the added 
feature of membrane bending coupling for shell problems; where the curvature 
is zero, there is no coupling between membrane, and bending effects, and TRSHL 
for such cases reduces to a combination element. The results for problems 
in this section have been obtained using TRSHL elements. 

Three buckling problems were investigated using the triangular plate and 
membrane elements. Following are the three different problems: (i) Buckling of 
a tapered column fixed at the base is shown in figure 3a. The area moment of 
inertia at any cross section can be expressed m the form 

CD 

where Ij is the moment of inertia at the top of the column (x = a) . Results 
for the buckling factor for a tapered column of I 1 /I 2 = 0.2 have been obtained 
from NASTRAN using TRIA2 and TRSHL, and an analytical solution from reference 
22 is given in table 1 for comparison, (ii) Buckling of a simply supported 
square plate subjected to uniform compression m one direction. Owing to 
symmetry, only one quarter of the plate (modeled with 2x2 mesh size) is used 
as shown in figure 3b. Results of the buckling factor from NASTRAN TRIA2 and 
TRSHL elements and the exact solution are shown in table 2. 

The nondimensional buckling factor X is represented by the formula: 



ir 2 P 

b 2 


( 4 ) 
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(m) The third problem considered is buckling of simply supported rectangular 
plate of aspect ratio a/b - 0.8 under m-plane bending loading shown m 
figure 3c. Due to symmetry, only half of the plate is used in the analysis. 
NASTRAN results using TRIA2 and TRSHL with different mesh sizes are shown 
m table 3, along with analytical results from reference 22. Table 3 clearly 
shows that the TRSHL elements gave a much better prediction of the critical 
buckling load than the TRIA2 elements. 
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191 




X 


(a) Tapered column. (b) Simply supported square 

plate under uniform 
compression . 




n ° r~ 


I 

1 

f " a 


1 

1 

b 






(c) Simply supported plate under 
in-plane bending. 


Figure 3. Column and plate geometry for TRSHL element buckling 
test problems. 
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Exact 


4.0000 














Table 3. Buckling factor for a simply supported rectangular plate 
of aspect ratio 0.8 under m-plane bending, v = 0.3. 
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APPENDIX E 

The NASTRAN source code subroutines that 
need to be modified or added for 
inclusion of TRIM6, TRPLT1 and TRSHL elements 
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APPENDIX E 


and 

The NASTRAN Subroutines that are modified to 
TRSHL elements are 

add the 

TRIM6, TRPLT1 

1 . 

DS1 

15. 

IFX7BD 

29. 

LD14 

2. 

DS1A 

16. 

LDjZil 

30. 

LD15 

3. 

EDTLZZZ 

17. 

LD^2 

31. 

LD21 

4. 

ELELBL 

18. 

LD03 

32. 

LD22 

5. 

EMGPR0 

19. 

IDM 

33. 

LD23 

6. 

GPTABD 

20. 

LD05 

34. 

LD34 

7. 

IFP 

21. 

LD06 

35. 

LINEL 

8. 

IFS1P 

22. 

LD07 

36. 

$FP1A 

9. 

IFX1BD 

23. 

LD*48 

37. 

0FP1BD 

10. 

IFX2BD 

24. 

LDjzS9 

38. 

tfFPSBD 

11. 

IFX3BD 

25. 

LD10 

39. 

0F1PBD 

12. 

IFX4BD 

26. 

LD11 

40. 

0F5PBD 

13. 

IFX5BD 

27. 

LD12 

41. 

SDR2B 

14. 

IFX6BD 

28. 

LD13 

42. 

SDR2E 


New Subroutines Added 

1. KTRM6S: Stiffness and mass matrix generation subroutine, single 

precision version, for element CTRIM6. 

2. KTRM6D: Stiffness and mass matrix generation subroutine, double 

precision version, for element CTRIM6. 

3. TL0DM6: Thermal load vector calculation for element CTRIM6. 

4. STRM61 Stress data recovery. Phase I for element CTRIM6. 

5. STRM62: Stress data recovery. Phase II for element CTRIM6. 

6. KTRPLS Stiffness and mass matrix calculations, single precision 
version (without the effects of transverse shear), for element CTRPLT1. 

7. KTRPLD Stiffness and mass matrix calculations, double precision 
version (without the effects of transverse shear), for element CTRPLT1. 
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8. TSPL1S. Transverse shear calculations, single precision version, 
for element CTRPLT1. This subroutine performs the numerical integration to 
obtain the contribution to the generalised stiffness matrix due to transverse 
shear effects. 

9. TSPL1D. Same as TSPL1S, double precision version. 

10. TSPL2S Calculations, single precision version, to obtain the matrix, 
[B 2 ] , relating curvatures to generalized coordinates (m the equation 

{x!> = [B 2 ] {a}). 

11. TSPL2D: Same as TSPL2S, double precision version. 

12. TSPL3S Calculations, single precision version, to obtain the matrix, 
[B i ] , relating transverse shear strains to the generalized coordinates (in the 
equation {yj} = [BjJ {a}) for use m the TSPL1S subroutine. 

13. TSPL3D* Calculations, double precision version, to obtain the matrix 
[Bi] , as in TSPL3S, for use m the TSPL1D subroutine. 

14. TL0DT1: Thermal load vector calculations m the absence of transverse 

shear effects for element CTRPLT1. 

15. TL0DT2: Numerically integrated contribution to the thermal load 

vector for element CTRPLT1 due to transverse shear effects. 

16. TL0DT3. Calculations to obtain the matrix, [Bj] * relating trans- 
verse shear strains to the generalized coordinates (m the equation 

(yi) = [Bi] {a}) for use m the TL0DT1 subroutine. 

17. STRPll: Stress data recovery, Phase I, for CTRPLT1 element. 

18. STRPTS: Calculations to evaluate matrices for recovery of shear 

forces for CTRPLT1 element. 

19. STRP12 Stress data recovery. Phase II, for CTRPLT1 element. 

20. KTSHLS Stiffness and mass matrix calculations, single precision 
version, for triangular shallow shell element CTRSHL. 

21. KTSHLD: Stiffness and mass matrix calculations, double precision 

version, for triangular shallow shell element CTRSHL. 

22. TL0DSL Thermal laod vector calculations for element CTRSHL. 
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23. STRSL1* Stress Data Recover/, Phase I, for element CTRSHL. 

24. STRSLV Calculations to evaluate matrices for recovery of shear 
forces for CTRSHL element. 

25. STRSL2- Stress Data Recover/, Phase XI, for element CTRSHL. 

26. DTSHLD . Differential Stiffness Matrix generation, double precision 
version, for element CTRSHL. 
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